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CLASSICAL FIFTH-, SIXTH-, SEVENTH-, AND 
EIGHTH-ORDER RUNGE-KUTTA FORMULAS 
WITH STEPS IZE CONTROL 


INTRODUCTION 

In two earlier papers [l ], [2], the author has described a Runge-Kutta 
procedure which provides a stepsize control by one or two additional 
evaluations of the differential equations. This earlier procedure, re- 
quiring an m-fold differentiation and a suitable transformation of the 
differential equations, yielded (m+4)-th order Runge-Kutta formulas - 
as well as (m+5)-th order formulas for the purpose of stepsize control. 
The stepsize control was based on a complete coverage of the leading 
local truncation error term. The procedure required altogether six 
evaluations of the differential equations, regardless of m. 

Since for m=0 no differentiations have to be performed, our earlier for- 
mulas represent, in this special case, classical Runge-Kutta formulas 
of the fourth order, requiring six evaluations of the differential equations 
and Including a complete coverage of the leading truncation error term 
for the purpose of stepsize control. 


In this paper we shall derive classical Runge-Kutta formulas of the fifth, 
sixth, seventh, and eighth order including a stepsize control procedure 
which is again based on a complete coverage of the leading local trunca- 
tion error term. Naturally, these new formulas require more evaluations 
per step of the differential equations than the known classical Runge-Kutta 
formulas without stepsize control procedure. 



However, they require fewer evaluations per step than the known classical 
formulas if Richardson's extrapolation to the limit is used for such for- 
mulas as a stepsize control device. Since the application of Richardson's 
extrapolation to the limit means practically a doubling of the computational 
effort for the benefit of the stepsize control only, it is worthwhile to look 
for a less expensive stepsize control procedure. 

Less expensive procedures have been suggested by different authors. 
However, these procedures generally do not make any effort to cover the 
truncation error, but rather try to estimate the truncation error from the 
last or the last few considered terms of the Taylor series. Since such a 
procedure has no mathematical base, these estimates are rather unreli- 
able. Generally, since the terms in a convergent Taylor series are de- 
creasing with increasing order, the last considered term will be larger 
than the leading truncation error term. Therefore, a stepsize control 
procedure based on the last or last few considered terms of the Taylor 
series will, in general, largely underestimate the permissible stepsize, 
thereby wasting computer time and building up unnecessarily large 
round-off errors. 


3. The new formulas of our paper contain one or more free parameters. By 
a proper choice of these parameters the leading term of the local trunca- 
tion error reduces substantially. This, in general, results in an increase 
of the permissible stepsize. This increase, together with the smaller 
number of evaluations per step, accounts for the superiority of our new 
formulas compared with the known Runge-Kutta formulas operated with 
Richardson's principle as stepsize control procedure. 

4. Naturally, the new classical Runge-Kutta formulas of this paper, being of 
the eighth or lower order, are in general less economical than our earlier 
Runge-Kutta transformation formulas [l], [2] which represent Runge- 
Kutta formulas of any desired order. However, our new formulas have 
certain advantages, since they require no preparatory work (like differen- 
tiation of the differential equations) by the programmer. The new formulas 
including the stepsize control procedure can easily be written as a sub- 
routine. 
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One might try to further raise the order of our new classical Runge-Kutta 
formulas hoping to make them still more economical. However, we be- 
lieve that we have more or less reached the optimum with our eighth-order 
formula. The examples that we ran show that the gain of our eighth-order 
formula RK8(9) if compared with our seventh-order formula RK7(8) is not 
very substantial any more. 

We made the same experience with a ninth- and a tenth-order Runge-Kutta 
formula that E.B. SHANKS has developed recently. These (not yet published) 
new formulas of SHANKS do not bring a substantial gain any more compared 
with SHANKS's eighth-order Runge-Kutta formula. 
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PART 1. FIFTH-ORDER FORMULAS 

SECTION I. THE EQUATIONS OF CONDITION FOR THE 
RUNGE-KUTTA COEFFICIENTS 


5. Although all formulas of this paper hold for systems of differential 
equations, we state them — for the sake of brevity — for a single differential 
equation 


y' = f(x,y) 


U) 


only. 


In the customary way, we set: 


f o = f ( x o> y 0 ) 

K-i 

f = f ( x o + «K h < y G + h Z 

X=0 


( K 



( 2 ) 


and we require: 


5 

y = y 0 + h Z c K f K + °< h8 ) 

fc— 0 


y = y Q + h Z V* + 0(h?) 

K=0 


(3) 


with h as integration step size. 

Equations (3) mean that we try to determine the coefficients a , p , 

K K A 

c c in such a way, that the first formula (3) represents a fifth -order and 

K K 

the second formula (3) a sixth -order Runge-Kutta formula. According to (2) 


4 



and ( 3 ) , this means that the coefficients a and P have to be the same in both 

K KX 

formulas for the first six evaluations of f. This restriction explains why we 
shall need eight evaluations (instead of seven) for our sixth-order formula. 

6 . Since the equations of condition for the Runge-Kutta coefficients, result- 
ing from Taylor expansions of (2) and (3) , are well known in the literature, we 
restrict ourselves to stating these equations. For Runge-Kutta formulas up to 
the eighth order these equations — in condensed form — are listed in a paper by 
J.C. BUTCHER ([3], Table 1). For the convenience of the reader, we list 
these equations for a sixth-order formula such as the second formula (3) in the 
customary summation form. However, since these equations then become some- 
what lengthy, we introduce the following abbreviation: 


_ X n X 
+ ^2“» + 


+ 0 ,a , 

KK - 1 K-l 


kX 


(k= 2,3 7; A = 1,2,3) 


(4) 


The 37 equations of condition for our sixth-order formula, listed in the same 
order as in BUTCHER'S paper, then read: 


TABLE I. EQUATIONS OF CONDITION FOR SIXTH-ORDER FORMULA 


(1,1) V c -1 = 0 

K = 0 K 

7 

(H, 1) Yea -7 = 0 
u . K K 2 

K-i 

7 

(hi, i) y c p -j = ° 

K K 1 6 


(m,2) ~ Yj 0 a 2 - ^ = 0 

2 LJ J K K 6 


2 i _ 


#C = 1 


dv.yij'jp 


2 K K 2 24 

K- 2 
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TABLE I. (Continued) 


^• 3) 2 * k “« p ki - i = 0 


k=2 


(IV - 4 > « 2 g X-i=° 


6 

(V,1) A 

k-4 

<v.2)t Z' s izV x p , 2 )- Ii5-° 

K — o \ A~ / 

<V,3> 

k=3 \ A-2 j 

<v ' 4) 5 2^>K3-lk = ° 

k-2 

< V ' 6) 22jV« P «2-» _# 

k=2 

<v - 7) I z’WiH 

k=2 

(V,8) | 

k=2 

< v - 9 > m Z' s ,“i-Ii5 =0 


K = i 


7c- 1 /X-i 

iMkv- 


1 

120 


(vi, i) E’UZ o 


k=5 


1 \rx A 

(VI ’ 2 > 5 2 \ 

k=4 


kA 




Jl 

x-i 


X=4 

Is ^L?2 


720 


= 0 


720 


= 0 
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TABLE I. 


(Continued) 


(VI, 3) 
(VI, 4) 
(VI, 5) 
(VI, 6) 
(VI, 7) 
(VI, 8) 
(VI, 9) 
(VI, 10) 
(VI, 11) 
(VI, 12) 
(VI, 13) 
(VI, 14) 
(VI, 15) 


E 

k=4 


^3 


1 ? T s Ci p p 

6 H L n K\ A3 / 720 
k=3 \a=2 / 


240 


= 0 


- 0 


E 

k=4 


A 

C 


K — 1 

E « 

A=3 


'A-i 


O' 


y i s p 

A/ , 


*a M .r'o A^ mi 


u=2 


180 


= 0 


a P „ 
kA A A2 


1 

180 
- = 0 


1 v 7 A / V 1 a 

2 Z/ c K l Z P 

k-3 \\=2 

5 J 3 C ,< (| 2^ Px1 )' I5 »° 0 

J_y 7 A p _J_ = 0 

24 Li * k-4 72ft 


24 44 k k4 720 
k=2 


y A 

L c 


K= 4 


Q! 

K K 


K- 1 /A-l 

E y( E 


1 y 7 a / 
O A C « 

2 K =3 ‘ *\ 


A=3 
K-l 

E ^x p 

A=2 


M= 2 
P 

kA A2 


am mIj 

/ 

1 


144 


7 . / 

'K-l 

S C a ( 

E * 

K=3 \ 

A=2 

1 y 7 a 



P 

6 Li k 

K k3 

k= 2 


7 

/K-l 

y c p i 

r v 

44 K K i 

i ^ 

k=3 ' 

( A=2 


144 

E ^x“x P x, ]-^ = 0 


kA A A1 
1 

— = 0 


— y c p p ~ — = o 

2 Li k k1 k2 72 ° 

K = 2 
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TABLE I. (Concluded) 


(VI, 16) 
(VI, 17) 
(VI, 18) 
(VI, 19) 
(VI, 20) 


— V c a 2 P - — = 0 

4 k k k k2 72 

k=2 

iy'c% p ! , - ~ = o 

2 k k k1 48 

K= 2 

« z' ? “ 3p , 

6 K K K 1 72 

k=2 

— t ! 

120 t-J 


K = 1 


Ac 1 

C Q! s - = 0 

(C K 720 


The Roman numerals in front of the equations in Table I indicate the order of 
the terms in the Taylor expansion. 

A similar table holds for a fifth-order formula such as the first formula 
(3) . We obtain this table from Table I by omitting the sixth-order equations 
(VI, 1) through (VI, 20), and replacing, in the remaining equations, c with c 

K K 

and the upper limit 7 of the K-sums with 5. 

Naturally, all equations of this new table for a fifth-order formula as 
well as those of Table I have to be satisfied simultaneously. 

7. The equations of Table I which represent necessary and sufficient con- 
ditions for a sixth-order Runge-Kutta formula, can, however, be replaced by 
a much simpler system of sufficient condition. 

First, we make the assumptions: 


a 5 = a 7 = 1 , cv 6 = 0; c t = Cj = 0, c 2 = c 2 , c 3 = c 3 , c 4 = c 4 , c 5 = 0, c 6 = £ 7 = c 5 (5) 
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which greatly facilitate the simultaneous solution of both systems of equations 
(for fifth-order and sixth-order Runge-Kutta formulas) . Then, we introduce 
the following assumptions (A) , (B) , (C) which are well-known to reduce our 
equations of condition to a great extent: 


(A) 


<B) 


(C) 


(*) 

P2i=f«i ( A2 > 

y*\ 3 . 2 _ 1 

(*) P51 ~ 2 as ~ 2 

(A5) 

(*) 

P 3 l = f«3 2 (A3) 

P 61 = ~2 = 0 

(A6) 

(*) 

P 41 = \ «! (A4) 

p _ J_ 2 = J_ 

(A7) 

<*) 

C2£21 + C3/3 31 + C4^ 41 + C 5 | 

!S; + ^h' u - ai,=o 

(Bl) 

(*) 

| 

C3£32 + c 4£42 + c 5 ' 
1 


( B2) 

(*) 

1 

c 4£43 + c 5 j 


(S3) 

( * ) 

C5 ' 

1 

\%+4 - ci<i - ai) 

(B4) 


c 5 ( £65 + £75 ) ~ c 5( 3 ” a 5) _ 0 

(B5) 



c 5£?6 “ c 5 

(B6) 

(*) 

C2 a 2£21 + C3«3£31 + C 4 04/3 41 + C 5 = 0 

(Cl) 


c 2 a 2£21 + C8 a i£31 + c 4 a 4£41 + c 5£71 = 0 

(C2) 


030:3/332/321 + 040:4(^42^21 + 043/331) + c 5(£?2/321 + £73/331 + £74/^41 
+ £ 75£ 5i + £ 7e£6i) = 0 • 


(C3) 



The asterisk (*) in front of some of trie equations indicates that these 
equations hold for the fifth-order formulas as well as for the sixth-order 
formulas. If the equations marked by an asterisk are split in two lines, the top 
line holds for the fifth-order formulas and the bottom line for the sixth-order 
formulas. The equations without an asterisk are required for the sixth-order 
formulas only. 

8. Inserting the assumptions (A) into Table I, one immediately finds the 
following identities: 


<m, 1) 3 (m, 2); (IV, 3) 3 3(IV, 4); (V, 7) 55 3(V, 9) ; 

(V, 8) = 6( V, 9) ; (VI, 8) = 2(VI, 7); (VI, 14) = (VI, 16); 
(VI , 15) = (VI, 17); (VI, 18) s 15(VI, 20); (VI, 19)= 10(VI, 20). 


( 6 ) 


Therefore, the equations on the left-hand sides of the identities (6) can be 
omitted from Table I. 

Using also assumption (B1 ) , three more identies are obtained: 

(IV, 1) 3 (IV, 2); ( V, 3) = 3(V, 4) ; < VI, 7) = 3( VI, 9) ; (7) 

thus eliminating equations (VI, 1) , (V,3), (VI, 7) from Table I* 

The assumptions (B) lead to the following identities: 

(IV, 2) = (111,2) - 3(IV,4); (V, 1) s (IV, 1) - (V, 5); \ 

(V , 2) = (TV ,2) - (V, 6) ; (V, 4) = (IV, 4) - 4(V, 9); I 

(VI, 1) = (V, 1) -(VI, 10); (VI, 2)= (V , 2) - (VI, 11); l (8) 

(VI, 4) = (V, 4) -(VI, 13); (VI,5)=(V,5) - 2 (VI, 16); l 

(VI, 6) = (V, 6) - 2(VI,17);(VI,9)3 (V, 9) -5(VI,20); I 

thus eliminating ten more equations from Table I. 
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Finally, we make use also of the assumptions (C) and obtain the following 
identities: 


(V, 5) s (V, 6); (VI, 3) = 3(V, 4) - 3(VI, 13); 
(VI, 10)= (VI, 11); (VI, 12) 3 3 (VI, 13); 

(VI, 16)= (VI, 17) . 


9. Cancelling all equations listed on the left-hand sides of the identities (6) , 
( 7) , ( 8) , ( 9) , Table I reduces to the following ten equations: 

(1,1), (11,1) (HI, 2), (IV, 4), ( V, 6) , ( V, 9) , (VI, 11), (VI, 13), 

(VI, 17), (VI, 20) . 



We arrange these remaining equations in two groups, (D) and (E) , the former one 
not containing any /3-coefficients. We obtain; 


(D) 


( * } |l °| +C2+C3 + C4+ jlJcsj = 1 

(*) c 2 a 2 + c 3 a 3 + c 4 a 4 + c 5 - | 

( * ) c 2 a j; + c 3 a| + c 4 a 4 + c 5 = ~ 

( * ) c 2 a 2 + c 3 a 3 + c 4 a| + c 5 = ^ 

( * ) c z a\ + c 3 af + c 4 a| + c 5 = ^ 
c 2 af + c 3 a 3 + c 4 af + c 5 = ^ 


(DO) 

(Di) 

(D2) 

(D3) 

(D4) 

(D5) 
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and: 


. P 52 1 

(*) c 2 a 2 P 22 + c 3 « 3 p 32 + c 4 a 4 P 42 + c 5 j p | - lg 


c 2 a 2 P 23 + C3Q-3P33 + c 4 a 4 P 43 + c 6 P 73 - 24 


C2«2 P 22 + c 3 a 3 P 32 + c 4 a 4 P 42 + C 5 P 72 “ jg 


C 3 a 3 3^ P S S + C 4 a 4 (0 42 P 22 + /?43 P 32> + 

c 5 (^72 P 2 2 + 3 ? 3 P 3 2 + ?7 4 P 42 + 37 5 P 52 + 37 6 P 6 2) = 


72 


Using (C) and introducing the abbreviations: 


P = p - B ffi 
f kA kX k 1 


(k = 2, 3, 4, 5, 6, 7; X = 1,2,3) 


(ID 


the last group of equations reduces to: 


(*) c 3 a 3 p 32 + c 4 a 4 p 42 + c 5 


P 52 / _ _L 


P 72 1 5 


(E) 


C 3 a 3 P 3 3 + c 4 a 4P43 + C 5P73 “ £4 


C 3 a 3 P32 + c 4 a 4P42 + C 5P72 ~ jg 


c 4 Q! 4^ 4 3P 32 + C 5(/ 3 73P32 + ^74?42 + 075P52 + 076P62) - ?2 


(El) 

(E2) 

(E3) 

(E4) 


Finally, the coefficients fl (k = 1, 2, 3, . . . , 7) are obtained from the 
well-known equations: 
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I 


(F) 


\ 


(*) /3 10 - a i (Fi) 


(*) fi 20 + @21 = a 2 

(*) @30 + @ 31 + @32 ~ a 3 

( * ) fiiO + @41 + @42 + @43 = “4 

(*) 050 + 051 + 052 + ^53 + 054 = a 5 

@60 + @61 + @62 + @63 + @64 + @65 

@10 + @11 + @12 + @13 + @14 + @15 


\ 



(F2) 


(F3) 


(F4) 

1 

(F5) 

o 

1! 

C£> 

B 

(F6) 

076 = a 7 = 1 

(F7) 


The asterisks ( * ) in ( D) , ( E) , ( F) have the same meaning as in equations 
(A) , (B) , (C) . In the following section we shall show how the sufficient 
equations of condition (D) and (E) can be solved, taking into account the 
assumptions (A) , (B), (C). 
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SECTION 1 1. A SOLUTION OF THE EQUATIONS OF CONDITION 
FOR THE RUNGE-KUTTA COEFFICIENTS 

10 . When solving our equations of condition, we shall try to express the 
8- and c- coefficients as functions of the o;-coefficients . However, not all four 
coefficients ax, as, a 3 , a 4 are independent of one another. In the following 
two numbers we shall show that a 4 and a 3 can be expressed by as . Therefore, 
we have only two independent a-coefficients; namely ai and as . Our ( 3 - and 
c- coefficients will then be expressible by ai and as only. 

11 . To find the relation between a 4 and as, we form: 


03(1 - a 3 ) P 3X + c 4 (l - a 4 )P 4 x . 


( 12 ) 


Using the assumptions (A), (B), (C), (D), we obtain from ( 12 ): 


c 3 (1 - a 3 ) 832 as + c 4 (1 - a 4 ) (0 4 2a 2 + 0 4 3 a 3 ) - — 


( 13 ) 


We further form the expression: 

as • (B 2 ) + at • (B 3 ) + a 4 -(B 4 ) - (El) , ( 14 ) 

thereby eliminating the term C5 (852 at + $53 at + 054 a 4 ) in (El) and obtaining: 

c 3 (1 - a 3 ) @32 al + c 4 (1 - a 4 ) (842 a2 + 04a <23) = gg • (10) 

We now eliminate the terms with c 3 in ( 13 ) and ( 15 ) and find: 

c 4 (l - a 4 ) 843 a 3 (a 3 - a 2 ) =-^- “ as . ( 16 ) 
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We establish a second relation between c 4 , / 3 43 and the a- coefficients. We 
start by forming: 


C 3 CI 3 P 3 I + C 4 «4 P41 + C5P51 

Using the assumptions (A), (B), (D), we find from (17): 

c 3 Qf 3/ 3 32 a 2 + C 4 CV 4 (/342Q - 2 + Pi 3 & 3 ) + C5i^52 a 2 + 053 a 3 + P 54 a 4) = ^ (18) 

By forming: 

Ci A t P3 1 + C4 0>4 P 4 1 + c 5 P51 ( 4 9) 

and using (A), (B), (D), we find in the same way: 

C3<*3 Ps2 a 2 + C 4 Q! 4 (/ 3 42 Q ! 2 + p i3 a 3 ) + c 5 (P 52 a 2 + p 53 a 3 + > 8540 : 4 ) = . (20) 

We finally form: 

(E3) -a 2 * (20) +a 2 • (18) - (El), ( 2 i) 

thereby eliminating the terms with c 3 and with c 5 . As second relation between 
c 4 , p i3 and the a-coefficients, we then obtain: 

1 J_ 

C4P43 a 3 Q 4 (03-0:2) (1 - <24) - 90 40 * ( 22 ) 


From the two relations (16) and (22) , we find as a relation between a 4 and a 2 : 



12. A further relation between a 2 , a 3 , and a 4 can be obtained from (Dl) , 
(D2), ( D3) , (D4) , and (D5) as condition of compatibility. Eliminating from 
these equations the terms with c 5 , c 2 , c 4 , c 3 , in this order, leads to: 


0?3 = 5 Q ?3 0/4 - 3(02 + 0/4) + 2 

10a2a;4 - 5 (q , 2 + 0/4)+ 3 


(24) 


Introducing (23) in (24) we find as a relation between a 3 and a 2 : 


Q!2 

“ 3 “ 15 as -lOas+2 


(25) 


13. The weight factors c (or c ) are now obtained from equations (Dl) 
(D2) , (D3) , (D4): K * 


A 1 

c 2 = c 2 


45a 2 - 30a 2 + 4 


180 (15a 2 - 10a 2 + 2) (5a 2 - 2) a 2 (a 3 - a 2 ) (“4 “ a 2 ) (1 - a 2 )\ 


a 1 

c 3 ~ c 3 _ 


15a 2 - 10 q! 2 + 2 


180 (5a 2 - 2) a 3 (a 2 - a 3 ) (a 4 - a 3 ) (1 - a 3 ) 


_ A 1 

c 4 ” C 4 


(5a 2 - 5a 2 + 1) ( 5a 2 - 2) 


(26) 


20 (15a 2 - 10a 2 + 2) a 4 (a 2 - a 4 ) (a 3 - a 4 ) (1 - a 4 ) 


A A 1 

c 5 = c 6 = c 7 = — - (c 2 a 2 + c 3 a 3 + c 4 a 4 ) 


14. From the two equations (El) , it follows: 


P?S =P 53 


(27) 
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Multiplying equations (B2), (B3), (B4) by a a at\, respectively, 
adding these three equations, and using (D2) and (D3), leads to: 

1 

C 3 P32 + C 4 P42 + C 5 Ps3 =— . (28) 


The three equations ( 28) , (El) , (E3) represent a system of three linear 
equations for the three unknowns P32, P42, P52. Its solution reads: 

t 3a 4 - 2 J 

P32 ~~ 180 c 3 (o !4 - a 3 ) (1 - a 3 ) I 

1 3q; 3 - 2 ^ 

P42 ~ 180 c 4 (cv 3 - a 4 ) (l-ar 4 ) / 

. 15o! 3 Q!4 - 12(0:3 + 0:4) + 10 1 

P52 = 180 ' c 5 (l - a 3 ) (1 - a 4 ) * I 


(29) 


15. The four equations (A2) , (A3), (A4), (A5) yield: 


a I 2 

P2i&l = 7«2 

1 2 

= Y«3 - P31 

1 9 

Pil a l=~£ a i ~ P41 


(30) 


The last equation ( 30) is obtained using equations ( B2) , ( B3) , ( B4) and 
equations (Di), (D 2 ). 
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Introducing (30) into equation (Cl) results in the following expression 

for p 41 : 


Pix — 


C 3 P31 ( 1 “ &3 ) ~ 24 


C 4 ( 1 - 014 ) 


(31) 


Since /3 32 is already known from the first equation (29) , equation (31) and the 
second equation (29) represent two equations for the two unknowns /3 42 , jS 43 . 
Their solution reads: 


c iPi2 ~ 


15a: 2 a 3 (o! 4 - ex 3 ) + 4(a 3 - a 2 ) - 6 q! 3 (q !4 - a 2 ) 
360a 2 (a 3 - a 2 ) (a 4 - a 3 ) (1 - a 4 ) 


i 5ct 2 - 2 

c 4^43 = " 720*0:3(0:3 - a 2 ) (1 - a 4 ) 


(32) 


From equations (B2), (B3), (B4), it now follows: 


c 5/^52 - c 2 (l - a 2 ) - C $22 - C4/3 42 \ 

c 5@53 = C 3 U - a 3) - c iPi3 \ (33) 

C5054 = C 4 (l - 0:4) 


It can easily be shown that the coefficients /3 52 , /? 53 , as obtained from 
(33), satisfy the last equation (29). 

Since the coefficients /3 41 a? 1 , /3 51 q: 1 are obtainable from 

(30) , all Runge-Kutta coefficients for our fifth order formula are now determined 
as functions of the parameters and a 2 - 
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16. We still have to find the coefficients P 6 ^and 0 7 ^ for our sixth- 
order formula which is required for the stepsize control. 

From the two equations (Ci) it follows: 


07i ~ 0si 


( 34 ) 


From (B6) we obtain: 


076 “ 1 • 


(35) 


Equation (B5) suggests: 


065 ~ 075 “ 0 . 


(36) 


From equations (A7), (El) , (E2) , the following expressions for p T1 , p^, p*, 
are obtained: 


C 5P71 


C 5P72 


C sP73 


1 . 60ag(l - g 3 )(1 - 0 : 4 ) + 2(304 -2) - 15q; 2 (1 - q^) 

360 a-z (1 - a 3 )(1 - < 24 ) 

1 . 15q?b - 12 ( a 3 +a4) + 10 

180 (l-« 3 )(l- 04 ) M 37 ) 


1 . 15(1 - ag )(1 - a?4)+ 15 a s a 3 a; 4 (1 - a 3 ) - 6a 3 a? 4 (l-o;g -a^ )- 4a;3 

360 (l-o' 3 )(l-0' 4 ) 
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Equations (37) yield the following solution for P 74 , 0 73 , P 7g : 


C5$74 ~ ~ en 


(12q a -5)(5 gj -5g a + l) 

60 (15a;i ~ lOafe +2)ct!4(a!2 - a 4 )(or 3 - a 4 ) 


05^73 ~ 


0's («3 


3 _ ®2) L 


60q?s q; 4 - 45q; 3 - 30q! 4 + 24 
360(1 -a*) 


c 5 074 0 ? 4 (0 ' 4 - Cte ) 


] 


\ 


(38) 


_ ft _ 1 f l5cfeQi4 - 12 (q! 3 + a 4 ) + 10 fl 2 R s i 

C5p7S "^T L 180(1 -a 3 )(l-« 4 ) - C b(P73«3-P74«4)J 


/ 


Finally, equations (B2), (B3), (B4) yield for the coefficients 
8 62 > 8 63 > 8s4 ■ 

c 5 862 = c 2 ( 1 “ CliS ) - C 3 P 3 3 - C 4 P 4 s ~ C5 P 7 2 
^5803 = C3 ( 1 - G! 3 ) - C 4 P 43 - C5 P 73 

c 5 P 64 = c 4 ( 1 - a 4 ) - c 5 P 74 


(39) 


17. It can be verified that the Runge-Kutta coefficients, as obtained in 
this section, satisfy indeed all our equations of condition (A) through (F), 
if c 0 , c 0 , 0i O , Pa, , , 0 7O are determined from (DO) and (F) . From the 
derivation of our Runge-Kutta coefficients it is not obvious that they 
satisfy equations (C3) and (E4), since these two equations have not been 
used in the derivation. However, when expressing all coefficients by 
as explicitly, it can be shown that the two equations (C3) and (E4) become 
identities in as . 
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SECTION III. THE LEADING TERM OF THE LOCAL 
TRUNCATION ERROR 


18. The leading truncation error term of our fifth-order Runge-Kutta formula 
is a sum of twenty members, each being of the form: 


T 


v 


[•••• 3v • h 3 


(40) 


The T are numerical constants that are characteristic for the individual 

v 

Runge-Kutta formula under consideration; they are, however, independent of the 
differential equation (1) . The bracket in (40) contains certain combinations of 
partial derivatives of the right-hand side of (1). The twenty T - values, 

expressed by the Runge-Kutta coefficients of Section II, can be obtained from 
the left-hand sides of equations (VI, 1) through (VI, 20) of Table I, if, in these 
equations, we replace c by c and the upper limit 7 of the k- sums by 5. 

K K 

If we list these T - values in the same order as in Table I, the first 

T read: v 

v 


T x = 0585484383283! 0?! 

Ts “ 2~ w K ^2) J ~ 


1 

720 

1 

" 240 


(41) 
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19. Naturally, it will be desirable to establish Runge-Kutta formulas with 

small T - values. 

v 

the sum of the terms (40) representing the leading term of the local truncation 
error can also be expected to be reasonably small. 


If we succeed in making all T - values sufficiently small, 

v 


Due to the special structure of our Runge-Kutta formulas (assump- 
tions (5)), the majority of our T -values are zero, as we shall see immediately. 


Because of (5), (27), (34), the equations of condition (C2), (D5), ( E3) 
hold for fifth -order Runge-Kutta formulas also. Since (D5) and (E3) correspond 
to equations (VI, 20) and (VI, 17) of Table I, we have T 20 = 0 and Tj 7 = 0. From 
T 20 = 0 we obtain, using ( 6) , ( 7) , ( 8) : T 19 = 0, T la =0, T 9 = 0, T 7 = 0, T g = 0. 
Correspondingly, from T 17 = Owe find, using (6), (8), (9): T 15 = 0, T e = 0, 

Tjg = 0 , t 5 = 0 , t 14 = 0. 

By a proper choice of the parameter Oj we can also make T X1 - 0. The 
straightforward computation, using the results of Section I and II, shows that 
Tjj = 0 requires 


2 5a 2 - 1 

a i = TiT ’ 


(42) 


From the identity (VI, 2) = (V, 2) - (VI, 11) in (8) , it follows that T X1 = 0 leads 
to T 2 = 0 also. Therefore, altogether, fourteen T -values of the leading trun- 
cation error term are zero. 

20. We now compute T^ using the results of Sections I and II. Expressing 
Ti 0 by the only free parameter as , the computation yields 

_ ! . (43) 
Tl ° 360 15 q| - 10a 2 + 2 
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and from the identity (VI, 1) = (V, 1) - (VI, 10) in (8), it follows: 

Ti = - T 10 . (44) 

21. We still have to determine T 3 , T 4 , Tis, Txs . Again, a straight- 
forward computation, using the results of Sections I and II, leads to: 

1 (3 cus - 1) (Seel - 5a!2 + 1) 

Tl3 = " 360 15 oj - 10a 2 + 2 ’ (45) 

Finally, from (8) and (9) we obtain: 


T4 T13 , T3 - 3Tj3 , T\s — 3 Ti 3 . 


( 46 ) 



SECTION IV. EXAMPLE FOR A FIFTH-ORDER 
RUNGE-KUTTA FORMULA 


22. From (43) through (46) , it follows that all T^-values become zero for 

oi 2 = 1/3. However, this value for a z leads to a% - 1, <*4 = 1, thereby making 
equations (Dl) through (D5) incompatible. Therefore, we have to discard the 
value a 2 = 1/3. This means that we cannot make all T^-values zero. This 

was to be expected, however, because otherwise we would obtain a sixth-order 
Runge-Kutta formula with only six evaluations. 

One might try to keep the T^-values small by choosing a 2 close to 1/3 

(this means cr 3 and a 4 close to 1) . However, such a choice of a 2 is not 
advisable, since this would lead to large values for some of the weight coeffi- 
cients (26) . 

23. Since we cannot make all T -values zero, we might try to make the 

v 

T -values (45) and (46) zero by assuming: 

v 

5«| - 5g!2 + 1 = 0 (47) 

This would leave us with only two T^ -values different from zero (Tj and T 10 ) . 
From (47) we would obtain: 


01 2 



(5 ±>fE) 


(48) 


and from (43), (44): 


T 10 


5 ±3>/5 
3600 


(49) 


Obviously, the minus sign in (48) would be preferable, since it leads to 
a smaller value Tjo . 
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However, the values (48) do not lead to a solution of our problem: for 
a z = (5 --/5) we would obtain from (47) : <23 = ~ (5 + -JH >) . From equations 

(Dl) through (D5) we derive, by eliminating c 5 , c 2 , and c 3 , the two relations: 


c 4 of 4 ( of 3 - a 4 ) (ot 2 - a 4 ) (1 - a 4 ) 

040:4(0:3 - Q! 4 ) (a 2 - 0:4) (i - a 4 ) 


= ~ [ 100 : 20:3 - 5(q: 2 + 0 : 3 ) + 3] 

= ^ [ 50 : 20:3 - 3(a 2 + «3) + 2 ] 


(50) 


For our values a 2 , ce 3 the right-hand sides of (50) would both become zero. 
Therefore, one of the following conditions must hold: 0:4 = 0, 0:4 = 1, 0:4 = a 2 , 
a 4 = 0:3 or c 4 = 0 . It can be shown, however, that for any of these conditions 
our equations of condition (A) through (E) lead to contradictions. 


24. However, by choosing for a 2 a value in the vicinity of ~ (5 - V~5) the 

terms (45) and (46) can be kept small. On the other hand, the terms (43) and 

(44) vary only slightly in the vicinity of as = yjj- (5 - J~5) for which value they 

have an extremum. Therefore, instead of choosing a 2 = ~ (5 - -/~ 5 ) « 0. 276 

4 

we might, for example, take a 2 = — « 0 . 267. This value is relatively close 

1 ^ 

to ~ (5 - V~5) ; it yields relatively simple Runge-Kutta coefficients — all weight 


factors (26) being non-negative — and the error coefficients T 10 and T 13 become 
relatively small. Table H shows the Runge-Kutta coefficients for a 2 = 4/15, 
as obtained from Section H. In this case, the leading term (the sixth-order 
term) of the truncation error obtained as the difference of the two formulas 
(3) , would read: 


TE _ 0 g (fo + f 5 - 4 ~ £ 7 ) h 


(51) 


From (43) and (45) we then obtain for our error coefficients: 


T ‘° = ' 2160 “ " 463 ‘ 10 " 3; T 


13 32400 


0.309*10 


-4 


(52) 
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TABLE H. RK5(6) 
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SECTION V. NUMERICAL COMPARISON WITH OTHER 
FIFTH-ORDER RUNGE-KUTTA FORMULAS 


25. For comparison, we consider some other fifth-order Runge-Kutta 
formulas. The first two fifth -order Runge-Kutta formulas ever developed 
are due to W. KUTTA ( [4] , p. 446/447) . Since both formulas of W. KUTTA 
contained minor arithmetical errors which later were corrected by this author 
( [5], p. 424), and by E. J. NYSTROM ( [6], p. 5), respectively, we state 
in Tables III and V KUTTA's formulas in their correct form. We also present 
in Tables IV and VI the error coefficients Tj through T 20 for KUTTA's formulas. 
These coefficients have been rounded to three significant digits. 


TABLE HI. KUTTA 1 
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TABLE IV. ERROR COEFFICIENTS FOR KUTTA 1 


n 

T 

T 

T 

T 

T 

5n+l 

5n+2 

5n+3 

5n+4 

5n+5 

0 

-0. 139- 10" 2 

-0. 556- 10 -3 

-0. 833- 10" 3 

-0.278- 10 -3 

0.278- 10“ 2 

1 

0. ill- 10~ 2 

0.417- 10" 3 

0. 833- 10" 3 

0. 139- 10 -3 

0. 139- 10- 2 

2 

0. 556- 10 -3 

0. 833- 10" 3 

0. 278- 10" 3 

-0.389- 10 -3 

-0.222- 10" 3 

3 

-0. 389- 10" 3 

-0. 222- 10- 3 

-0. 250- 10" 3 

-0. 167- 10“ 3 

-0. 167- 10" 4 


TABLE V. KUTTA 2 
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TABLE VI. ERROR COEFFICIENTS FOR KUTTA 2 


n 

T 

5n+l 

T 

5n+2 

T 

5n+3 

T 

5n+4 

T 

5n+5 

0 

-0. 139- 1(T 2 

0 

-0. 833- 10 -3 

-0.278- 10 -3 

0. 278- 10 -2 

1 

0 

0.417- 10 -3 

0. 833- 10" 3 

0. 139- 10“ 3 

0. 139- 10" 2 

2 

0 

0. 833- 10" 3 

0. 278- 10 -3 

-0. 556- 10" 3 

-0. 926- 10 -4 

3 

-0. 555- 10" 3 

-0. 926- 10 -4 

-0. 278- i0 -3 

-0. 185- 10" 3 

-0. 185- 10" 4 


Comparing Tables IV and VI with the error coefficients (52) of our 
Runge-Kutta formula, we notice that in both formulas of KUTTA the coefficient 
Tjo is three times as large as in our formula, and the coefficient T 13 even 
nine times as large as our coefficient T 13 . Furthermore, in KUTTA's formulas 
the term T 5 , which is zero in our formula, is six times as large as our 
largest term T t = -T 10 .Finally, KUTTA's first formula has twenty non- 
vanishing error coefficients T , and his second formula has still seventeen 
v 

such terms, whereas our formula carries only six non-vanishing error 

coefficients T . 

v 

For these reasons we may expect that our formula (Table II), having 
a smaller truncation error, shall integrate a problem in fewer steps than 
KUTTA's formulas (Table III and Table V) . 

26. More recently, several other authors (H. P. KONEN, H. A. LUTHER, 

J. A. ZONNEVELD, etc. ) have derived more fifth-order Runge-Kutta formulas. 
However, we cannot consider all these fifth-order formulas, especially since 
most of them are really not much of an improvement with respect to the error 
coefficients if compared with KUTTA's formulas. 

27. For the numerical comparison of our formula RK5(6) with KUTTA's 
formulas, we apply these formulas to an example of a system of two differential 
equations, for which the exact solution is known. Therefore, in our example, 
we can easily check on the accumulated errors of these formulas. 
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Example: 


y ' = - 2xyTogz, z'=2xz-logy 

Initial values: x = 0, y = e, z 
o J o 

„ , , cos(x 2 ) 

Exact solution: y = e , z 


= 1 
o 

sin(x 2 ) 
= e 


( 53 ) 


Since no other step-size control procedure seems to be known for KUTTA's 
formulas, we have applied RICHARDSON'S extrapolation to the limit as 
stepsize control for these formulas. 


All computations were executed on an IBM-7094 computer (16 decimal 
places) . The same tolerance (10 -i6 ) for the local truncation error was allowed 
for all formulas. 


Table VH shows the results of our computer runs for example (53); all 
runs cover the interval from x = 0 to x = 5. 


TABLE Vn. COMPARISON OF FIFTH -ORDER METHODS FOR EXAMPLE (53) 


- 


Results for x = 5 and Tolerance 10“^ i 

A* ! 









Number of 
Substitutions 
per Step 

Number 
of Steps 

Total Number 

Running Time 
on IBM- 7094 

Accumulated Errors in y and z 

Method 

of Evaluations 

(min) 

Ay 

Az 

KUTTA i 

11 

6335 

G9685 

7.10 

-0.3071. 10” 12 

-0.2772 . io-1 2 

KUTTA 2 

11 

6290 

69190 

6.95 

-0.2383- 10“ 12 

-0,2802. 10" 12 

RK5(6) 

8 

4779 

38232 

3.84 

+ 0.1072* 10~ 12 

-0,2190. 10" 12 


The table shows that in this example our formula RK5(6) requires approximately 
55 percent of the evaluations or 55 percent of the computer running 
time compared with KUTTA's formulas, all formulas yielding about the same 
accuracy. 

Such a saving could be expected since the computational effort per step — 
including stepsize control — is smaller for our method than for KUTTA's 
formulas. Furthermore, we save steps since our local truncation error is 
smaller. 
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PART 1 1. S IXTH-ORDER FORMULAS 

SECTION VI. THE EQUATIONS OF CONDITION FOR THE 
RUNGE-KUTTA COEFFICIENTS 


28. For sixth-order Runge-Kutta formulas with stepsize control, we 
allow for ten evaluations of the differential equations per step. 

Therefore, we have instead of (2) and (3) : 


f Q = f ( x o> y 0 ) 


k-1 


f = f(x + a h, y + h V 0 f ) (k = 1,2,3, ... ,9) 
K 0 K O LJ- kA A 

A— u 


(54) 


and 


y = y + h V cf + 0( h 7 ) 

J J 0 K K 

K - 0 


A 

y 


y + h Yj c f + 0(h 8 ) 


K = 0 


K K 


(55) 


29. Since Table I contains the Taylor terms up to the sixth-order only, we 
need a supplementary table for the seventh-order terms. Such a table can again 
be obtained from J. C. BUTCHER'S paper ([3] , Table I). It would contain 48 
equations. Because of the length of this supplementary table, we do not reproduce 
it here, but refer directly to BUTCHER'S paper. We denote the equations of this 
new table — in the same order as in BUTCHER'S table — by (VII, 1) through (VH.48). 

This supplementary table in combination with Table I contains all 
equations of condition for a seventh-order Runge-Kutta formula. Naturally, 
in Table I the K-sums now run up to 9 (instead of 7) . Since we are Rooking 
for a pair of Runge-Kutta formulas (55) , the original Table I (with c replaced 

by c ) must also be satisfied for our sixth-order formula. 
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30. We proceed in a very similar way as in Part I and make the following 
assumptions: 


a 7 = a s = 1, oe 8 = 0, c 7 = c 7 = 0, c 2 = c 2 = 0, c 3 = c 3> c 4 - c 4 , c 5 - c 5 ,i 

A A „ A A f 

Cg = Cg, c 7 = 0, Cg - c 9 - c 7 


@31 “ Pil ~ ^51 ~ @61 - ^71 - @ 81 - ^91 - 0 
Further, analogous to Part I, the following assumptions shall hold: 


(*) P 2i = f «2 

(Al- 2 ) 

(*)P 61 =fa| 

(A 1-6) 

(*) P 31 

(Al-3) 


(A 1-7 ) 

(*)P 4 1=|«4 2 

(Ai-4) 

*81 =£«» = <> 

(A 1-8) 

(*) 1*51 “ 2 0 i 

(Al-5) 

1 2 1 
p 91 - 2 a 9 - 2 

(Al-9) 

(*) p 32 = 3 

(A2-3) 

. . , 1 3 1 

(*)P« = 3 a l = g 

(A2-7) 

(*) P 42=f a 4 

( A2-4) 

P 8 2 = £ «8 = 0 

(A 2 - 8 ) 

(*) P 52 = 3 “5 

(A2-5) 

i a 1 

p 92 = 3«9 = 3 

(A2-9) 


(*> P 62 = T a 6 (A2-6) 
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(B) 



(*) c 2021 + c 3031 + c 4041 + c 5051 + c 6061 + c 7 | ^ + ^ 

(*) + C4/342 + C5/352 + Cg^g 2 +c 7 I ^ + 


(*) 
(*) 
(*) 
( * ) 


04043 + C5/353 + Cg^g3 + C7 


j^73 i 
/ /^sa + £93 5 


+ °» | pJJ + fi M | 

c Jfw I 

|086 + 096 j 


C7 (087 + 097) 


c 7 098 


Cl(l 

- a t ) 

= 0 

(Bl) 

C 2 (i 

- a 2 ) 

= 0 

(B2) 

C 3 (l 

- a 3 ) 


(B3) 

C 4 (i 

- a*) 


(B4) 

c 5 (l 

- " 5 ) 


(B5) 

C 6 U 

- «e) 


(B6) 

c T (l 

- a 7 ) 

= 0 

(B7) 

Cj(l 

-Q! 8 ) = 

C7 

(B8) 


(C) 


(*) 


c 3^3032 + c 4 Q! 4042 + C 5 Q! 5052 + c e a 6062 + c 7 jjg^j = ^ 

C3Q! 3/332 + 040(4^42 + C5a|/3 5 2 + C0Of0/302 + C7 3g2 = 0 

c 4 “4043032 + c 5 a 5 ( 053032 + 054042) + c eP t 6(063032 + 064042 + 065052) 

+ 07 (^ 93/332 + fi 94042 + 095052 + 096 062 + ^97^72 + 0 98082) = 0 


(Cl) 
( C2) 


(C3) 


The asterisk ( * ) in front of some of these equations means now, natur- 
ally, that these equations must hold for the seventh-order formula as well as 
for the sixth-order formula. 
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31. We now reduce the necessary and sufficient equations of condition of 
our supplementary table to fewer and simpler sufficient equations of condition by 
inserting the assumptions (Ai), (A2), (B),and (C) into the table. 

Since the reduction is completely analogous to the procedure of Number 
8, we shall present here the results of the reduction only, without going into 
all details. The reduction leads to equations (VII, 24), (VII, 29), (VII, 37), 
and (VII, 48) as the only four independent equations of the table. All other 
equations of the table can be expressed by these four equations and by some 
sixth-order equations of Table I. 

For the benefit of the truncation error terms, which we shall consider 
later, we now list those equations which can be expressed by (VII, 24) and 
(VH.29): 


and 


(VII, 1) s (VI, 1) -(VII, 21) 
(VII, 2) 3 (VI, 2) -(VII, 22) 
(VII, 3) = (VI, 3) - 3 (VII, 24) 

(VII, 4) 3 (VI, 4) - (VII, 24) 

(Vn,21)3 (VD.24) 

(VH, 22) 3 (VII, 24) 

(VH,23) = 3 (VII, 24) 


(VH,5) = (VI, 5) 
(VII, 6 ) 3 (VI, 6 ) 
2(vn,7) 3 (VI, 8) 
(VH, 8) 3 (VI, 8) 

(VH, 9) 3 (VI, 9) 
(VH,25) 3 4(VH,29) 
(VII, 26) 3 4 (VII, 29) 
(VH,27) 3 3(vn,29) 
(VH,28) 3 6(vn,29) 


4 (VII, 29) 
4 (VII, 29) 
6(VII, 29) 
6 (VII, 29) 
(VII, 29 


(57) 


(58) 


The remaining equations of our supplementary table can be expressed 
by equations (VII, 37) and (VII, 48). However, because of our assumptions (56), 
these two equations do not contribute to the truncation error of our sixth-order 

formula . 
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It should be noticed that the fifth and sixth identities in ( 57 ) are correct 
only if assumption (C 3 ) is assumed to hold. Since condition (C 3 ) does not hold 
in the case of a sixth-order formula, it can not be followed from (VII, 21 ) - (VII, 24 ) 
and (VII, 22 ) = (VII, 24 ) that the truncation error coefficients Tax and Tss are 
identical with Ta4 . 


32 . We now list for our pair of sixth- and seventh-order Runge-Kutta formulas 
the sufficient equations of condition that result after we have performed all the 
reductions in Table I and in the supplementary table. 

We arrange them again in two groups (D) and (E) in the same way as in 
Part I. 



and 


A° > + C3 + C4 + C5 + Cg + 

O ' 

C30 3 + c 4 o 4 + c 5 Q! 5 + c 6 o 6 + c 7 = ■£ 

C3CU3 + c 4 a| + C50 5 + c 6 a| + c 7 = ^ 

03(23 + 04014 + c 5 a| + CgOg + c 7 =~ 

C3O3 + C4O4 + c 5 a| + CgOg + c 7 = -g 

C3O3 + 0404 + C 5 05 + CgOg + = ~ 

C3O3 + 0404 + C5O5 + CgOg + C7 = — 



(DO) 

(Dl) 

(D 2 ) 

(D 3 ) 

(D 4 ) 

(D 5 ) 

(D6) 
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(El) 


(*) C3Q! 3P33 + c 4 a 4 P43 + C5QI5P53 + c 6 a 6 P 63 + c 7 


? 73l __L 

? 93( 24 


C 3 ® 3 P 34 + c 4 a 4 P 44 + c 5 a 5 p^ + c 6 a 6 p 64 + c 7 P M = — 


^ ^ \ C 3 Q; 3 P 33 + C 4 a 4 P 43 + C 5 C 25 P 53 + CgagPg 3 + c 7 P 93 - 2Q 


1 4^4/^431^33 + C5CV 5 (/353P33 + /354P43) + CgQ!g(^ 03 P 33 + ^g 4 P 4 3 + /305P53) 

1 


+ C 7 (/I93P33 + /I94P43 + 35-^53 + ^96 P 63 + ^97 P 73 + 


140 


(E2) 
( E3) 


(E4) 


Group D is an obvious extension of the former group (D) of Number 9, 
the last equation (D6) representing equation (VII, 48) of the supplementary 
table. In group (E) the first equation represents equation (E2) of the former 
group (E) of Number 9. Equations (E2), (E3), (E4) of group (E) correspond 
to equations (VII,29), (VII, 37), (VR, 24) of the supplementary table. 


Because of our assumptions (A2) the equations (El), (E3) of our former 
group (E) in Number 9 become identical to the equations (D4), (D5) and can, 
therefore, be omitted from our new group (E). Introducing our new assumptions 
(A2) and (Cl) into equation (E4) of our former group (E), it follows immediately 
that this equation becomes identical with equation (El) of our new group (E). 

Since the extension of our former group ( F) for sixth-order Runge- 
Kutta formulas is quite obvious, we omit this extension here. 
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SECTION VII. A SOLUTION OF THE EQUATIONS OF CONDITION 
FOR THE RUNGE-KUTTA COEFFICIENTS 


33. In the case of our sixth-order Runge-Kutta formula, six a -coefficients, 
a^, a 2 , CH3, 014, a 5 , a 6 , are at our disposal. We try again to express the 
/3- and c-coefficients as functions of these a -coefficients. Again, we shall 
show that not all six a -coefficients are independent of one another. Obviously, 
from (Al-3) and (A2-3) , it follows 


«2 



(59) 


34. A further relation between the a -coefficients is obtained from equations 
(Dl) through (D6) as condition of compatibility. Eliminating from these 
equations the terms with c 7 , c 6 , c 3 , c 4 , c 5 , in this order, leads to: 


«5 


i 

7 


35a 3 a 4 a 6 - 21(a 3 a 4 + a 3 Q! 6 + ae 4 a 6 ) + 14(a 3 + a 4 + a 6 ) - 10 
10a 3 a 4 a 6 - 5(a 3 a 4 + a 3 a 6 + a 4 a 6 ) + 3(a 3 + a 4 + a 6 ) - 2 


(60) 


35. We now establish a relation between a 3 , a 4 , oiq , which will 
allow us to express a 2i a z , a 5 by a 4 and a e so that a 4 , a 4 , a 6 will be the 
only free a-coefficients. To find such a relation, we proceed in the 
following manner: 

From (Al) we find by multiplying these equations with c 3 (l - a 3 ) , 
c 4 (i - a 4 ) adding them and using (D2) and (D3): 


c 3 ( 1 - ct 3 ) P 31 + c 4 (l - a 4 ) P 41 + c 5 (l - a 5 ) P 5i + c G (l - a 6 ) P 6 i = ^ (61) 

and in the same way from (A2) , using (D3) and (D4) : 

c 3 ( 1 - Q! 3 ) P 32 + C 4 (l - a 4 ) P 42 + c 5 (l - a 5 ) P 52 + c 6 (l - a 6 ) P e2 = ~ (62) 
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A third equation is found by forming: 

c^(B2) + ajj(B 3 ) + a|(B 4 ) + ae|(B 5 ) + ar^(B6) - (El) , ( 63 ) 

thereby eliminating in ( 63 ) the terms with c 7 . We obtain: 

£3(1 - a 3) P33 + 04(1 - a A ) P43 + c 5 (l - a 5 ) P53 + c 6 (l - a? e ) P 6 3 = ( fi 4 ) 

From equations ( 61 ) , ( 62 ) , ( 64 ) , we can easily eliminate the terms with c 3 
and C4. We obtain: 

Cg ( 1 - a 5 ) /?54« 4 («4 - as) («4 - a 2) + c 6 (l - a 6 ) [P Gi a l( a i - a 3) (a 4 - a 2 )j 

( 65 ) 

( 66 ) 

( 67 ) 

( 68 ) 

and find: 


+ ^ 65 ^ 5 (a 5 - a 3 ) (a 5 - a 2 ) ] = [5a2a 3 - 2(a z + a 3 ) + 1 ] 

If we had multiplied (Al) with c 3 (l - o§) , c 4 (l - a\) , ... , we would have 
obtained instead of ( 61 ): 

c 3 (l -a§) P 31 + c 4 (l -a\) P 41 + c 5 (l -0:5) P 51 +c 6 (l - af) P 61 = 
and by multiplying (A 2 ) with c 3 (l - at 3) , 04(1 - a|) , ...» 

c 3 (i - a 3 ) P32 4 c 4 (l - &l) P 42 + c 5 (l - <*5) Pg2 + c 6 (l - afg) ^62 = 3g 
Analogous to ( 63 ) we now form: 

a|(B 2 ) + « 1 (B 3 ) + al(B 4 ) + ajj(B 5 ) + <*g(B6) - (E 3 ) 
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C 3 (l - «i) P 33 + C 4 (l - o 4 ) P43 + c 5 (l - 05) P53 + c 6 (l - o|) P 63 


70 


(G9) 


From equations ( 66 ), (67), (69), we can again eliminate the terms with 
c 3 and c 4 and find: 


c 5 (l - 0*5) £5404(04 - o 3 ) (o 4 - o 2 ) + c 6 (l - al) [^ 64 Q' 4 (a 4 - o 3 ) (o 4 - o 2 ) 
+ P 65 a 5 (a 5 - a 3) (a 5 - a 2 ) ] = [84020:3 - 35 (o 2 + o 3 ) + 18 ] . 


( 70 ) 


Obviously, equations (65) and (70) permit the elimination of the terms with 
c 6 and yield: 


c 5^54( a! 6 “ <*5) U “ a 5 ) ^4(^4 " ® 3) (®4 “ ® 2) 


1 

2520 


{21(1 + a 6 )[5a20 3 


- 2 (a 2 + a 3 ) + 1] - 2 [84 o 2 o 3 - 35(a 2 + a 3 ) + 18]} . 


(71) 


To eliminate C 5 / 354 , we need a second relation between 05/354 and the o -coefficients. 
We start from (Al- 5 ) and (A2-5). Multiplying these equations with c 5 and 
eliminating the terms with /3 53 yields: 


c 5^54 a 4( Q! 4 “ ®3) “ c 5^52® 2 (® 3 ~ a 2) 


1 ? 

- c 5 05(3o 3 - 2o 5 ) . 


( 72 ) 


From- (B2) , (Cl), (C 2 ) we find by eliminating the terms with /3 e2 and /3 72 : 


c 5/ 3 52( Q! 6 - or 5 ) ( 1 - 05) = -C3/3 32 (o 6 - O3) (1 - 03) - C4/342 (o 6 - O4) (1 - O4) (73) 


Combining (72) and (73) yields: 
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c 5^54 ( ^6 " a 5) (1 " a 5) a 4( a 4 " a 3) CaP 32 a 2( a 3 " “2) ( a 6 “ a 3 )(*■ a a) 

- c^^2 a 2( a 3 - a 2 ) (“6 - “ 4 ) (1 - «4)- gC 5 a|(3a 3 - 2a 5 ) (a! 6 -Q! 5 ) (l-a 5 ) 


( 74 ) 


From (Al-4) and (A2-4) we find by eliminating the terms with $43: 

1 ? 

/ 342 « 2( Q! 3 “ “2) = “ a|( 3 a 3 - 2a 4 ) 


(75) 


and from (Al-3) and (59) 


^Z2 a 2( a 3 ~ a 2) 


= -|a^(3a 3 - 2a 3 ) 


(76) 


Inserting (75) and (76) into (74), we obtain: 


05/354 (a 6 - a 5 )(i - a 5 ) a 4 (a 4 - a 3 ) - - j a 3 [c 3 af(a 6 - a 3 ) (1 - a 3 ) 


+ c 4 a 4 (a 6 - a 4 ) (1 - a 4 ) + c 5 a|(a 6 - a 5 ) (1 - a 5 )] 

+ ^ [c 3 a|(a 6 - a 3 ) (i - a 3 ) + c 4 a|(a 6 - a 4 ) (1 - a 4 ) 


+ c 5 aij(a 6 - a 5 ) (i - a 5 )] 


(77) 


The values of the brackets in (77) can easily be obtained from equations (Dl) 
thr ough (D6) by eliminating c 7 and c 6 . For the first bracket we obtain 

^•(5a e - 3) and for the second one ^(3a e - 2) . Inserting these values in 

(77) , we find: 
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05/354 (a 6 - a 5 ) (i - a 5 ) a 4 (a 4 - a 3 ) - [ 2 ( 3 a 6 - 2 )- 3 a 3 ( 5 a 6 - 3 )] 


(78) 



From equations (71) and (78) we can now easily find the wanted relation between 
the a -coefficients by forming: 


(78) • (a 4 - a 2 ) - (71) = 0 (79) 

We obtain: 

42a 4 a 6 - 28a 4 - 21 q!6 + 15 

a 3 =~ (80) 

15a 40:6 - 9au - 6a@ + 4 


36. Having subjected the a -coefficients to the restrictive conditions (59) , 
(60) , and (80) , we can now easily find the c- and /J -coefficients in terms of 
the a -coefficients . The procedure is quite similar to that of Part I. 


From equations (Dl) through (D 6 ) we obtain for the c-coefficients: 


A 1 . 

C 3 = C 3 = iT 


10a 4 a 5 a 6 - 5(a 4 a 5 + a 4 a e + a 5 a 6 ) + 3(a 4 + a 5 + a 6 ) - 2 
< 23 ( 0-4 - a 3 ) (a 5 - a 3 ) (a 6 - a 3 ) (1 - a 3 ) 


A 1 
C4 — C4 — 


60 


10a 3 a 5 a 6 - 5(a 3 a 5 + a 3 a 6 + a 5 a 6 ) + 3(a 3 + a 5 + a 6 ) - 2 
a 4 (a 3 - a 4 ) (a 5 - a 4 ) (a 6 - a 4 ) (1 - a 4 ) 


a 1 

c 5 = C 5 = 


60 


10 a 3 a 4 a 6 - 5 (a 3 a 4 + a 3 a 6 + a 4 a 6 ) + 3 (a 3 + a 4 + a 6 ) - 2 
a 5 (a 3 - a 5 ) (a 4 - a 5 ) (a 6 - a 5 ) (1 - a 5 ) 


(81) 


A 1 

Cg — Cg - 


60 


10a 3 a 4 a 5 - 5 (a 3 a 4 + a 3 a 5 + a 4 a 5 ) + 3(a 3 + a 4 + a 5 ) - 2 
< 2 6 (<23 - a 6 ) (“4 - a 6 ) (®s - “6) “ Q e) 


A A 1 

c 7 = c 8 = c 9 = - - (c 3 a 3 + c 4 a 4 + c 5 a 5 + c 6 a 6 ) 
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37. 


From (Al-2) and (59), it follows: 


3 21 Qfx = ' 


2 

«3 


and from (Al-3) and (59): 



From (Al-4) and (A2-4) , we find: 


1 al(3a 3 - 204) 
342 = "e" * q- 2 - q/2 ) 

1 aH 3ois ~ 2 a*) 

04B -~g ' Q/g (q; 3 - q/ 3 ) 


( 82 ) 


(83) 


(84) 


Equations (B2), (Cl), (C2) can now be considered as three equations for 
052» j3g2» /5 72 , since /3 32 and p i2 are already known. From these equations, we 
obtain: 


c 3 (a 6 - 0:3) (1 - o; 3 ) / 3 32 + c 4 (q! 6 - a 4 ) (1 - or 4 ) / 3 42 


052 “ 


c 5 (ce 6 - or 6 ) (1 - a B ) 


c 3 (Q !5 - a 3 ) (1 - a 3 ) 0 32 + c 4 (a 5 - a 4 ) (1 - a 4 ) 0 42 


062 ~ 


c 6 (a 5 - a 6 ) (1 - ar 6 ) 


(85) 


072 ~ ~ ~~ ( c 3032 + c 4042 + c 5052 + c 6^62) 
c 7 
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From equations (Al-5) and (A2-5) , we now obtain j3 53 and £54: 


P 53 ~ 


P 54 


1 ■> 

- ag(3a 4 - 2 a 5) - p 52 a 2 (a 4 - a 2 ) 
a 3 (a 4 - a 3 ) 

a^(3a 3 - 2a 5 ) - j 3 52 a 2 (a 3 - a 2 ) 
a 4 (a 3 - a 4 ) 


1 2 


( 86 ) 


We eliminate P 73 from equations (El) and (E3) and find: 


Pe3 c 6 a 6 (l - a 6 ) 


— - c 3 a 3 (l - a 3 ) P 33 - c 4 a 4 (l - a 4 ) P 43 


- c 5 a 5 (l - a 5 ) 



(87) 


Since P 33 , P 43 , P 53 are already known from (83), (84), (85), (86), P e3 is a 
Jknown quantity now. 


Equations (Al-6) , (A2-6) and the definition of P 63 can now be considered 
as three equations for jS 63 , /3 65 , since /3 62 is already known. Their solution 

reads: 


P 6 ia 4 a 5 - P 62 (a 4 + <*5) + p e3 “ / 3 62 a 2( Q! 4 - a 2 ) (a 5 - a 2 ) 


. P 63 “ 


a 3 (a 4 - a 3 ) (a 5 - a 3 ) 


“ P 62( Q! 3 + ^5) + P 63 “ PG2 a z( a 3 ~ a z) ( a 5 “ ^2) 


Pii ~ 


a 4 (a 3 - a 4 ) (a 5 - a 4 ) 


( 88 ) 


^l«3 a 4 " ^2(«3 + a 4) + 1^3 “ / 3 62 Q! 2(Q ! 3 " «2> («4 " «2> 


Pg5 ~ 


o' 5 (a 3 - a 5 ) (a 4 - a 5 ) 
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the quantities P ei , P 62 , Pg 3 being represented by the right-hand sides of 
equations (Al-6), (A2-6), (87). 

From equations (B3) through (B6) , we obtain the coefficients /3 73 , ^74, 
/3 75 , /3 76 in the form: 


c 7 /3 73 = c 3 (l - a 3 ) - {cJ3 i3 + C5/3 53 + Cg/3g 3 ) 


c 7 / 3 74 = c 4 (i - a 4 ) - (05^54 + Cg/ 3 64 ) 


c 7^75 = 05(1-05) - Cg/ 3 g 5 


C7/376 - Cg ( 1 - Og) . 


( 89 ) 


38. We still have to determine the coefficients j3^ and /3^, required for 
our seventh-order formula. 

We determine first the coefficients fl . From (B8) , it follows: 

yA, 

P30 = 1 » ( 9 °) 


while (B7) suggests: 


8g? — 837 — 9 • 


(91) 


From the two equations (Cl) , we obtain: 


892 — 872 , 


(92) 


and from the two equations (El) : 
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The still missing coefficients /? 93 , 0 94 , /S 95 , 0 96 can be obtained from 
equations (Al-9) , (A2-9) , the definition of P 93 , and the definition of P 94 . We 
find: 


P 9l“ 4 “5“6 - P 32 (a 4 g 5 + a 4 a 6 + asdj) + P 93 (a 4 + a 6 + a 6 ) - Pm -/39 2 a 2 (a 4 - (“s ~ a 2> < a 6 a 2> 

093 = a 3 (Q ' 4 - ex 3 ) ( oi 5 - <x 3 ) (a 6 - a 3 ) 

PjjasQ’sO'e - p 92( a 3“5 + «3 Q 6 + + P 93< a 3 + »5 + a 6> " P 94 " 092 a 2( g 3 ~ a 2> ~ a 2> < tt 6 ~ g 2> 

094 = Ci((a3 - ff 4 ) («s - «<) (“6 - a «) 

PajagO^ag ■ p 92( a 3 a 4 + “3 a 6 + a 4 a G) + ^93 ( a 3 + «4 + a 6 ) - % - ffa2tt 2 ( g 3 ~ ( a 4 ~ a 2) ( a 6 ~ a 2> 

095 = _ < 25 ( 1*3 - a 5 ) («4 - “ 5 ) ( g 6 - “ 5 ) 

P 91 g 3 g 4 a 5 " Ik( a 3 a 4 + “3 a 5 + a 4 a 5> + ^3( a 3 + «4 + a 5> " % ' 09Z“2( a 3 “ a l> ( a 4 ~ O' 2) ( a 5 ~ g 2) 

096 = g 6 ( a 3 * g 6) (“4 - “6)(“5 ' a 6) 


(94) 


The quantity % is obtainable from (E2) as: 


P94 - ~ ( - c 3 a 3 P 34 - c 4 o; 4 P 44 - c 5 Q! 5 P 54 - c e cv 6 P 64 ^ 


(95) 


The coefficients j3 can now be found from equations (B) : 
8A 


082 “ 072 “ 092 ~ 0 


083 - 073 “ 093 


084 - 074 “ 094 


085 “ 075 “ 095 


086 _ 076 " 096 


087 - 9 


(96) 


All Runge-Kutta coefficients for our pair of sixth- and seventh-order formulas 
are now known. Under the restrictive conditions (59), (60), and (80) our 
coefficients satisfy all equations of condition (Al) , (A 2 ) , (B) , (C) , (D) , and 
(E). 



SECTI ON V 1 1 1. THE LEAD I NG TERM OF THE LOCAL 
TRUNCATION ERROR 


39. According to our supplementary table we have 48 terms T (t> - 1,2,3,. . 48) 

that contribute to the leading truncation error term of a sixth-order Runge-Kutta 
formula. We obtain these terms from the left-hand sides of the equations of our 
table by replacing c^ with c^ and the upper limit 9 of the K-sums by 7. For- 
tunately, because of our assumptions (56) and our conditions (Al) , (A2) , (B) , 
and (C) , the majority of the terms are zero. For those T^, that are not 

zero, we obtain from (57) and (58): 


Tj = - T 21 , T 2 — T 22 , T3 — 3T 2 4, T4 — T 2 4, T 2 3 - 3T24 


(97) 


T 5 = - 4T 29 , T e — 4T 29 , T 7 — 3T 29 , T 8 — 6T 29 , T 9 — T 29 
T 2 s = 4T 29 , T 28 = 4T 29 , T 27 = 3T 29 , T 28 = 6T 29 


(98) 


with: 


l 


„ , , c a 

6 K K 

K = 5 


1 V 


7 

I 

k=4 
7 

I 

K = 4 


K - 1 / A-i 

~K ~ 1 / A-l 


840 


840 


K-l 


T 24=7 51 c a [ T, P 

" 6 u . k kI / j „ kA ) 

K = 4 \ A=3 


1 


A3 840 


(99) 


_1_ y 1 

T 29 24 \ °K a K P K4: 840 

k-3 
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40. The expressions (99) can be simplified by further use of the equations 
of condition. We finally obtain: 


T 29 - 24 c 7( r V4 - p 94) 


T 24 “ g c 7^83^33 + ^84 P 43 + /% P 53 + ^86-^63) f 

\ ( 100 ) 

T 21 - T24 - ~ C 7 (/3 83 /3 32 + /3 8 4/?42 + /3s&^52 + ^86/ 3 62) a 2 l 

3 a l I 

T 22 - 9 “Z ( T 24 ~ T 2l) + T 21 . I 

Z « 2 / 


We observe that only T 22 depends on aj. By a proper choice of o'! we might 
be able to keep T 22 small. 
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SECTION IX. EXAMPLE FOR A SIXTH-ORDER 
RUNGE-KUTTA FORMULA 


41 . It would present an extremely laborious and unpleasant problem to express 
all c- and /3-coefficients by the free coefficients a j, a 4 , a; 6 and to determine 
these a -coefficients in such a way that the error terms (100) become reasonably 
small. 


However, since we have to consider these error terms to obtain an effective 
sixth-order Runge-Kutta formula, we used an electronic computer to determine the 
error terms T 29 , T 24 , T 21 , T 22 / a l for a variety of pairs of coefficients a 4 
ot 6 . 

We selected ct 4 and a 6 from the following set: 1/2, 1/3, 2/3, 1/4, 

3/4, using all positive fractions with denominators 2, 3, 4, 5, 6, 7, 

8, 9, 10, 11, 12, 14, 15. Inspecting the print-out of the above-mentioned error 
terms, we choose for our Runge-Kutta formula a combination of a 4 and a e that 
leads to reasonably small error terms. 

We choose: 

i. 6 1 

= ~2 ’ ( a t = 2 a d- ( 101 ) 


This combination leads to: 


T 29 « +0. 376- 10“ 5 , Tjy « -0. 919- 10" 4 , T 21 « -0. 273- 10 -4 , 
T 22 w -0. 758- lO" 4 


( 102 ) 


The other non-zero error terms of our formula can then be computed from 
(97) and (98) . 

The Runge-Kutta coefficients of our formula were determined from the 
equations of Section VII and are listed in Table VHI. 
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TABLE Vm. RK6 (7) 



00 . 


o 

tH 


CO 


LO 

CO 


Truncation Error Term: TE = + £? -f 8 -ig )h 


(103) 
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SECTION X. NUMERICAL COMPARISON WITH OTHER SIXTH-ORDER 

RUNGE-KUTTA FORMULAS 


42. The first sixth-order Runge-Kutta formulas were developed by A. HU'Ta 

I 7] . [8]. We use his second formula ([8] , p. 23), since it has somewhat simpler 
coefficients. HUTA's formulas are based on eight evaluations of the differential 
equations per step. 

Later, J. C. BUTCHER [9] and other authors found sixth-order Runge- 
Kutta formulas that require seven evaluations per step only. We consider one 
of BUTCHER'S formulas ([9] , p. 193) as of special interest, since this formula has 
favorably small truncation error coefficients. 

We computed the error coefficients T^ (u = 1,2,3, . . . ,48) for HUTA's 

formula as well as for BUTCHER'S formula. Comparing them with our error 
coefficients (102) , we found that HUTA's largest coefficient T^ is more than 
32 times as large as our largest coefficients T 3 = - 3T24 and T 23 = 3T24, while 
BUTCHER'S largest coefficients T 6 and T 26 are less than three times as large 
as our coefficients T 3 and T 23 . While our formula RK6(7) of Table VIH has 
only 18 non-vanishing error coefficients T^, HUTA's formula has 36 such 

coefficients and in BUTCHER'S formula all 48 coefficients are different from 
zero. 

43. For the numerical comparison of our formula RK6(7) with HUTA's and 
BUTCHER'S formulas, we again apply these formulas to our problem (53) . We 
ran our formulas on the computer in exactly the same way and under exactly 
the same conditions as in Number 27, again using RICHARDSON'S principle as 
stepsize control procedure for HUTA's and BUTCHER'S formulas. Table IX 
shows our results. 

In this example our formula RK6 (7) requires only approximately 58 
percent of the evaluations and of the computer running time compared with 
BUTCHER'S formula. If we compare our formula with the earlier HUTA 
formula, this percentage decreases even more to approximately 40-42 percent. 
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TABLE IX. COMPARISON OF SIXTH-ORDER METHODS FOR EXAMPLE (53) 
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PART III. SEVENTH-ORDER FORMULAS 

SECTION XI. THE EQUATIONS OF CONDITION FOR THE 
RUNGE-KUTTA COEFFICIENTS 


44. For a seventh-order Runge-Kutta formula RK7(8j with stepsize control we 
allow for thirteen evaluations of the differential equations per step: 


h = f(Xo, y 0 ) 

K-l 

f K = f(x -Fo^h, y D + h E P a f x ) ( K = 1,2,3, 

X=o 

10 

y = y 0 +h ck£k + ) 

12 

y = y Q + h E c K f K + 0(h 9 ) 

K = 0 


12 ) 


(104) 


(105) 


45. We now need the equations of condition for the eighth-order terms. These 
are also listed in BUTCHER'S paper ([3], Table 1). There are 115 such 
equations, and we shall refer to them as equations (VTII, 1) through (VIH, 115) - 
in the same order as in BUTCHER'S paper. 


46. To reduce these 115 necessary and sufficient conditions to a system of 
simpler and fewer sufficient conditions, we make - very similar as in Part I 
and in Part II - the following assumptions: 


Qtio = a is = 1, an =0; c x = c x = 0, c 3 = 03 = 0, c 3 = c 3 = 0, 

c 4 = c 4 = 0 , c 5 =c B , Q? = q, , o, =0,, c% = c 3 , c 9 = &, , 

C10 “ C U - C12 - c 10 

P31 = 04 i = 3 5i = Pei = P71 = Psi = Pai = P101 ' Pm = Pia = 0 

Pss = P32 = P?2 = Ps 2 = Pas - Pios = P112 = Pias = 0 

and: 


\ 

>-( 106 ) 


/ 
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(Al) 

1 8 

Pvi 2 

(v = 2,3,4, . 

12 ) 

(Al-v) 

(A2) 

P V 8 =T “V 3 

(v =2,3,4, • 

12 ) 

(A 2 -v) 

(A3) 

1 4 

( v = 5 , 6 , 7, . 

12 ) 

(A3-v) 


XI C V$V3 + C 10 

V=& 

| Pl03 \ 

\ Pll3 + Pl23 ' 

= c 3 ( 1 - CK 3 ) = 0 

(B3) 


9 

X ° V ^V 4 + C 10 

V=5 

| 3 104 i 

V 3 H4 + 3l34 ' 

= C 4 (l-G!4) = 0 

(B4) 


9 

X c v + Cl ° 

V = 3 

( 3 106 \ 

v Pi 15 + Pl26 

= Cs(l-G! 6 ) 

(B5) 

(B) « 

9 

X C V^V 6 + C 10 

V— 7 

/ 3 106 \> 

v 3 116 + 3lSB ' 

= Cs (l-«e) 

(B6) 


9 

C y 0v7 + Cio 

/ 3 107 l 

\ 3ll7 + 3l37 f 

= C 7 (l-tv 7 ) 

(B7) 


Cg &gg + ClO 

/3ios l 

l 3ll0 + Pl28 J 

= Cg ( 1 “ Ota ) 

(B 8 ) 


Cio 

/ 3 ice \ 

l 3ll9 + 3l39 / 

= C9 (l-CKg) 

(B9) 


Cio 

( 3 1110 + 3l31o) 

= C 10 (l-Q!lo) = 0 

(BIO) 


^ C 10 

3 1211 

= Cio (1-Qfll) = C 10 

(Bll) 
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r 


9 


(Cl) 


E 

LjL = S 


C p“p P p 4 +C l° 8 124 = 0 


(C2) 


(C) -< 


E C jJL. ^(JL. 8 |J,3 I” C 10 0123 — 0 

|i = 5 


P= 


u*=e 


v 


E 

P=8 


(C3) 


Pjjpi + CioPi 24- 

= 0 


(C4) 

V 1 \ 
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+ c 10 

X! Pl2v3y3 = 0 

(C5) 

‘ V=6 ' 




'M- 1 V 


11 



+ Cio 

X] Pl2V0\H = 0 

(C6) 

V=8 / 


v=£ 



47. The reduction leads to equations (VIII, 57), (VIII, 68), (VIII, 86), and 
(VHI, 115) as the only four independent equations for the 8-th order terms. 
After the reduction we again arrange the remaining independent equations 
of condition in two groups (D) and (E) in the same way as before: 


(D) 


«.♦{£} - * 


10 


V=5 


C\j Oi 


v ~ \L + 1 


(|i — 1,2,3, 7) 


(DO) 

(Dp) 
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_1_ 

35 


(El) 





9 /y-1 \ 1L Pl °V P V3 

X c p Op,( X Pp,vPv3 j +C 10< V=i 

H=s \v=4 / 1 II 

[ X ^ i2v ^ v3 

' V=4 


X! C p, ^p, Pp- 5 + C W Pk 5 


X C P-^ S PpL4 +C loPl3i = 


3|J,V FV 4 ) + C X o XI 3 X2 V P\)4 

/ V=5 


1 

140 


_1_ 

240 


(E2) 

(E3) 

(E4) 

(E5) 


Obviously, equation (D7) represents equation (VIII, 115) of BUTCHER'S table, 
while equations (E3), (E4), (E5) are identical with BUTCHER'S equations 
(VIII, 68), (Vni, 86), and (VHI,57). Equations (El) and (E2) originate from 
seventn-order terms and, under our new assumptions, are identical with 
equations (E2) and (E4) of Part II. Equations (El) and (E3) of Part II may 
be omitted now, since, because of our new assumptions (A3), they now 
become identical with (D5) and (D6). 
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SECTION XI I. A SOLUTION OF THE EQUATIONS OF 
CONDITION FOR THE RUNGE-KUTTA COEFFICIENTS 


48. From equations (Dl) through (D7) we obtain a relation between a 5 , o% , 
dr, , a B , Og as condition of compatibility. The relation can be written in the 
form: 


= 70q% otg Q% Q ?9 -42(0% a 2 Q%+. . .+o; 7 9 k )+28(q % q; 7 +. . .+ 0 % q 9 )-2 0(q%+. . . +ofe)+15 

2 70a% a : 7 o% 0 % -35(0% o; ? g%+. . .+ 0*7 a% a% )+21(a% o; 7 +. . .+ 0 % o% )-14(o%+. . .+Gi 9)+10 


This rather lengthy expression simplifies to at = 1/2, if we require a% + a% =1 
and a s + a 9 = 1. Therefore, we assume for the following: 

<*5 ae + a 7 = 1, o; 8 + ofe = 1, 

since such a choice of a% , o% , a 7 , a 8 , Qfe also leads to rather simple ex- 
pressions for the weight coefficients c 5 , Cs , . . . , c l0 . We obtain: 


16 14(54-50% +1)4 -14(54 -5o% +1) o% +(144 " 14o% + 3)^ 
= 105 ' (l-2o%) 3 (l-2a 8 ) s 


1 7a e (l-o:a )-l 

420 o%(l— c% )(l-2a%) s (a%-G%)(l-a%-0! 8 ) 

1 7 q%(1-q; 6 )-1 

420 ’ ot B (1-0*8 )(l-2o:8 ) s (cife -o% )(l-a 6 -Qte ) 


1 70o% (l-a%)o* 8 (1-Q% )-7[q* 6 (1-o%)+o% (l-a%)] + l 

420 o; 6 (1-0% )a% (1-0% ) 


> 


J 


Naturally, we have assumed in (109) that o% , a 7 , o% , o% are different from 
one another and different from 0, 1/2 , 1. 


(107) 


(108) 


(109) 
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49. Further restrictive conditions between the ^-coefficients shall now be 
established. 

From (A 1-2) and (A2-2), it follows: 

2 

Qfi =“ oi s (110) 


and from (Al-3) and (A2-3): 

2 

OLZ =~ a 3 ( 111 ) 


From (A1 5), (A2 5), and (A3-5), we obtain the following relation between 
cfe and 0 : 4 : 


or because of a 5 = 1 / 2 : 


1 4a 4 - 3 as 

“ 3 “ 2 “ 5 3 ^ T - 2«5 


1_ # 80:4 - 3 

8 3of 4 - 1 


( 112 ) 


( 112 a) 


From (108), (110), (111), (112a), it follows that all a -coefficients can be 
expressed by 04 , 0 % , a Q . 

Later in No. 52, however, we shall show that a B is already determined 
by our equations of condition and by our assumptions. Therefore, a 4 and 
ofe shall be the only free ce-coefficients. 


50. We now determine the B -coefficients. From (Al-2) we obtain: 


and from (Al-3) and (111): 






(113) 


(114) 
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Equations (A 1-4) and (A2-4) lead to: 

2 

1 CV4 3oi3 -20!4 
042 - 6 a 3 -as ’ as 

1 a '4 3 as -2o?4 

P43 - - 6 Cl 3 -013 ’ a 


(115) 


and equations (Al-5) and (A2-5) lead to: 


$53 — ’ 


2 

0'S 

a±-a 3 


P54 =- — 


2 


g a4;-0>3 


3 0?4 — 2 O 5 
a>3 

3 0:3 -2 as 

« 4 


(116) 


Equations (A 1-6), (A2-6) and (A3-6) represent three linear equations for the 
three unknowns 3 e 3 , 3 6 4 , and 3 65 . Their solution reads: 


2 

l 6o4 0fB -4(a4 t Q;5)ofe + 3afe 

063 " 12 " 6 a 3 (oi 4 ,- a 3 )(as~ o> 3 ) 

1 60:3 0's -4(q! 3 + 0' 5 )cK 6 +3ofe 

064 = 12 "°^ oujo3 -~a^)(a 3 - «4 ) 

12 

]_ Gc^ G ?4 - 4(0 ! 3 + <¥ 4)05 +3ai6 

065 ~ 12 o ?5 (g ?4 -a s )(a 3 - as) 


\ 


> (H7) 


J 


Assuming: 


3?3 - 0 


(118) 


we can determine, in the same way, 374 , £ 75 , Ptb from equations (Al-7), 
(A2-7), (A3-7): 


2 

„ 1 2 6ff5Qte -4(0' 5 +0fe)O7 +3a!7 

“74 - 12 0(7 aj _ (0'5 - al ) (cf 6 -"a! )" 

1 6q'4 0' s -4(q:4 +Q' 6 )CV7 +3q?7 

075 = 12 017 Gfe" (aT - O eTR c* -as) 

1 a 6 q!4o;5 -4(04 +q; b)®7 + 3qv 

076 = 12 ^ ^(~a4~~a 6 T(^’5~ :: '^T _ 


\ 

>■ (H9) 

/ 
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From equations (B3), (Cl) and (C3), we can now determine Pea, 3® , and 
^103 = 0123 : 


083 

093 = “ 

0103 = 0123 


Cs(Gfe - Q?5)(l ~®5)053 +c g (Qjg ~ Q! 6 )(1“06 )0©3 
McUs-OaKl-O's) 

C 5 (as - 0' 5 )(l-ff 5 )P S 3 + C 6 (a a - Ofe)(l-Q! 6 )P63 
C3 (o 8 -0' 9 )(l-O' 9 ) 

C 5 (q;8-Q , 5)(C«s -a 5 )053 + C G (dfe-Q! 6 )(Q!9-0; 6 )%3 


A 


> ( 120 ) 


c l£> (1 — Og )(1“0!9 ) 


Assuming: 


0104 — 0124 , 


( 121 ) 


we obtain, in the same way, P 84 , 094 , and Pi^ = Piai from equations (B4), (C2), 
and (C4): 

n c s (q?9 -a 5 )(l-Q' 5 ) 054 + Ce (Qfe -0te)(l-06)Po4+0 ? (09 -Q? 7 ) (1-0(7)374 ^ 

084 


09 4 


C 3 ( 0:9 - 0 's )( 1 - Qfe ) 

c 5 (ae -Qi5)(l— Q' 5 ) 054 +c 6 (a 3 - oie )(1~«6 )06 4 +Q? («a -Qfr )(l-o? 7 )0?4 
Cq (o: 8 -a 9 )( 1 -Qfe ) 


>( 122 ) 


05 (<% -a 5 )(a9 -«5 )054 + Ce (Cte - 0 's )(Qfe -a 6 ) 3 6 4 +C 7 (Oe -«7 )(CH9-0'7 )074| 
Pl04: _ Pl24 


c 10 (1 - Qfe )(l - 0 ; 9 ) 


J 


51. From equations (B) we obtain by multiplying equation (B3) with at, 
equation (B4) with at, etc. , adding all these equations and using equations 
(D4) and (D5): 

C5 P54 + Cg Pg4 + C7 P74 + C 3 P 0 4 + Qa P94 + Coo Pl 04 = gQ ( 123 ) 

This equation together with equations (El) and (E4) can be considered as a 
system of three equations for the three quantities Pa 4 , P 94 , P 104 , since 
P 54 , P&4 , P 7 4 are already known from No. 50. 

The system has the following solution: 
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c a («9 -QteHl-Qte ^ Ps4 = 840 _3 )" c ^ (°b _Qr s HI-otbJPs* -c, 3 (oj-at)(l-a 6 )P G 4 -c? (cfc - 0*7 )(1-07)P 7 4 

00(09 -a Q )(l-09) P j4 = ( 4 og ~ 3 )-c 5 (o 8 -os )(l-o 5 )P & i -c^ (Oc-o G )(l-o e )Pg4 -07(00 -o^)(l-o 7 )P 7 4 

Ci 0 (1-Qfe )(l*Ob )Pio4 [28oa o>9 -24(oa4 Qfe )+2l]-c 5 (ofe-Os )(Oa -05 )Ps4 -c 3 (ofe-o 6 )(oa -Qfe )P 64 

-O7 (0^“0*7)(0 8 -<*7)P 7 4 ^ 


( (124) 


The quantities P^ , P 64 , P 74 which occur on the right-hand sides of equations 
(124) can be expressed by a* and a& . Using the results of No. 49 and No. 50 
we find: 


P-= 4 -' 3 - 3 h [8 “ J < 9 “» 8 - 10 ^ 


Pt4 = 24“ (l-Qte) a C(6c^-5Qfe-4a 6 + 3) + 0:4 (140^-100^ + 2)] 


J 


52. Equations (Al-8), (A2-8), (A3-8) and the definition of P^ may be written 
in the form: 


Ps5 + Ps6 Ofe + Pgf7 Otr/ 
3 85 tts + Pae Qfe + Per? Gb? 


Pcs ®5 + Pe6 + Pet? 0?7 


Pes Gt% + Pee Cfe + Pgr? 0?7 


= 2 Qfe - Pea Q!3 - Pa4 CK 4 

= “g* - Pea - Pa4 °S, 

— 0?8 - P83 Q§ - Ps4 

~ P84 - Pe3 0?3 - P&4 OI 4 




>(126) 

J 


Since we know 833, 6 a4 , P34 from equations (120), (122), (124), we can con- 
sider equations (126) as a system of four linear equations for the three 
unknowns 8 3S , 89 s , && . 

As condition of compatibility, we obtain from (126) the following relation: 
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(127) 


12 


ae [6o(5Qte <*> - 4a 0 (a 5 Qfe + as a? + Qfe a 7 ) + 3cv| (a^* Qfe + otn )] -Pe4 


- Pa3 a 3 («5 —Qfa )(Qfe -of 3 )(o *7 -ce 3 ) -8 e4 a 4 (ces -fv 4 )(a e -a^)(a 7 -a * ) = 0 


In (127) we now have to insert our expressions for P 84 , 8® , and P 84 and express 
all 3- and c -coefficients by the cv-coefficients. Finally, we have to insert our 
restrictive conditions (108) and(112a)for the a -coefficients, thereby converting 
(127) into a condition between o% and cts since the terms with a! 4 drop out. 

A tedious but straightforward computation finally reduces the left-hand side of 
(127) to a product of several factors. One factor is a function of a B only and can 
be made zero by a proper choice of a B . The other factors contain a 6 also. 

They can not be made zero for rational values a G , a e , if we assume that a 6 , a 7 , 
a s , % be different from another and different from 0, 1/2, 1, in order to pre- 
vent the denominators of (109) from becoming zero. 

If we omit all factors that are different from zero, we are finally left with: 

6a B - 7« s + 2 - 0 . (128) 

From (128) follows ae = 1/2 or afe = 2/3. Since 0's = l/2has again to be dis- 
carded, we obtain: 



(129) 


as resulting from condition (127). 

The restrictive conditions (103), (110), (111), (112), (129) reduce the independent 
a -coefficients to two, namely cc* and Oq . 


53. Since the computation of the remaining 8-coefficients and auxiliary ex- 
pressions is straightforward, we just indicate from which equations they are 
obtained without stating the explicit expressions for them: 

Pas, Pae > Pa? from (A 1-8), (A2-8), (A3-8) 

P 94 , Pio* = P 124 (El), (E4) 

P95, P96 > P 97, Pss from (A 1-9), (A2-9), (A3-9), P 94 
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0106 , Pice ) 01(77 j 0io8 > 0ice from (B5), (B6), (B7), (B8), (B9) 

0125,01%, 01%, 0i®, 0i®. from (A 1-12), (A2-12), (A3-12), (C5), (C6) 

0i2io = 0, 0i2ii = 1 from (BIO), (Bll) 

0113 = 0, 0H4 = 0 , 0115, 0116 ) from ( (B3), (B4), (B5), (B6) 

0117 , 011B , 0119 , 01110 = 0 j ( (B7), (B8), (B9), (BIO) 

All Runge-Kutta coefficients for our pair of seventh- and eighth-order formulas are 
now known if we still determine c 0 , c 0 from equation (DO) and the coefficients 
0 vo ( v = 1 , 2 , 3 , . . . , 1 2 ) f : rom the extended equations (F ) . If we submit the 
a-coefficients to the restrictive conditions (108), (110), (111), (H2), (129), 
the Runge-Kutta coefficients satisfy all our equations of condition. 
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SECTION XIII. THE LEAD ING TERM OF THE LOCAL 
TRUNCATION ERROR 


54. There are 115 terms T v ( v = 1, 2, 3, . . . , 115) that contribute to the 
leading truncation error term of a seventh-order Runge-Kutta formula. However, 
because of our assumptions and our conditions (Al), (A2), (A3), (B), and (C), 
most terms Tv are zero. Indeed, there are only 40 error coefficients Tv 
that are different from zero. All those 40 coefficients T v are multiples of 
four coefficients T^g , T52 , Tg? , and Tgg : 



55. The four basic error coefficients T« , Ts, Tg? , T® can be written in the 
following form: 


10 

/ K— 1 

X-1 

" UL-X 

/V-l v-1 

T -® = E Ck “k: ( 

E 

E 

E 

E PvpPpx) 


\x = 5 

H=4 

_ V=3 

\p = 3 /J 




Tea 


1 

= 120 ^ CkQ?k Pks " 


1 

5760 


1 

5760 
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SECTION XIV. EXAMPLE FOR A SEVENTH-ORDER 
RUNGE-KUTTA FORMULA 


56. For our seventh-order Runge -Kutta formula we select for the independent 
parameters < 24 , 0,3 such a pair of values for which the four error terms (131) 
become small. 

Since it seems hardly possible to find such a pair of values on, a e analytically, 
we again ran a variety of combinations a> 4 , Qfe on an electronic computer and 
printed out the error terms (131) for each combination. Inspecting these 
printouts, we chose for our Runge-Kutta formula the following combination: 

= 5/12, ofe = 5/6 . (132) 

This combination results in rather small values for our error terms (131): 

T «~ 0.164«10“ 5 , TssW -0.567* 10" 6 , Tg? w 0. 765* lO -6 , « -0. 919- 10' (133) 

It also leads to relatively simple Runge-Kutta coefficients, which are 
listed in Table X. 
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41 

Truncation Error Term: TE = (f 0 + fio - fu - fis )h 


(134) 
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SECTION XV. NUMERICAL COMPARISON WITH OTHER 
SEVENTH-ORDER RUNGE-KUTTA FORMULAS 


57. The only other 7-th order Runge-Kutta formula, known in the literature, 
is due to E.B. SHANKS ([lO], p.34). It is based on nine evaluations of the 
differential equations per step. 

Computing the error coefficients T v (v = 1,2,3,..., 115) for SHANKS's formula, 
we found that these coefficients are rather large compared with the coefficients 
(130), (133), of our formula RK7(8). 

While our formula RK7(8) contains only 40 non -zero error coefficients T v , 
SHANKS's 7-th order formula has 102 non-zero terms T v . The largest co- 
efficient T v in SHANKS's formula is -0.353* 10 -3 , and there are 15 terms T v 
in SHANKS's formula which are larger in absolute value than lO” 3 . The 
largest coefficient T v of our formula RK7(8) is -0.459* 10 -5 and in our formula 
there are 16 coefficients T v which are larger in absolute value than 10 -6 and 24 
such coefficients which are smaller in absolute value than 10 -8 . 


58. For the numerical comparison of our formula RK7(8) with SHANKS's 
formula we applied these formulas again to our problem (53) in exactly the 
same way as in Part I and Part II. We again used RICHARDSON'S principle 
as stepsize control procedure for SHANKS's formula since no other satis- 
factory stepsize control procedure seems to be known for SHANKS's formulas. 
Table XI shows the results of our comparison. 

TABLET XI. COMPARISON OF SEVENTH-ORDER METHODS FOR EXAMPLE (53) 




Results for x = 5 and Tolerance 10 _ls 

A. . 

Method 

Number of 

Number of 
Steps 

Total Number 

Running Time 

Accumulated Errors 

Substitutions 

of Evaluations 

on IBM-7094 

in y and z 



per Step 

(min) 


Az 


SHANKS 

17 

1423 

24 191 

2,49 

-0. 1332* 10 _1? 

-0.7377-10' 13 






-la 

-la 

RK7(8) 

13 

818 

10 634 

1.12 

-0.2509.10 

-0.5135* 10 
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The table shows that, in this example, the computer running time for our 
formula RK7(8) is only about 45% of the running time of SHANKS's 7-th order 
formula. The relatively large number of integration steps (1423), required 
by SHANKS's formula, reflects the unfavorable magnitude of the error co- 
efficients of SHANKS's formula. 
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PART IV. EIGHTH-ORDER FORMULAS 

SECTION XVI. THE EQUATIONS OF CONDITION FOR 
THE RUNGE-KUTTA COEFFICIENTS 


59. An eighth-ordsr Runge-Kutta formula RK8(9) with stepsize control can 
be established if we allow for seventeen evaluations of the differential 
equations per step: 


fo = f(Xo , Yo ) 

K-l 

f K = f(xo + a; K h, y 0 + hI^0K\fx) 


14 

y =y o +hZc K f K + 0(h 9 ) 
16 

y = yo + hScKf«: + o(h 10 ) 

K* 0 


(135) 

(K = l,2,3,..., 16) I 

(138) 


60. Since BUTCHER'S paper ([3], Table I) lists the equations of condition 
up to the eighth-order only, we have to find the ninth-order equations that we need 
for our stepsize control by either extending BUTCHER'S table or by using a 
scheme suggested by SHANKS [10 ] . 


We like to mention, however, that more recently another approach, based on 
quadrature formulas and also applicable to Runge-Kutta formulas of any order, 
was described by C. STIMBERG ([11], pp.1/9). For fourth-order Runge-Kutta 
formulas, C. STIMBERG has presented his method in more detail in another 
paper [12], For a similar approach based on quadrature formulas see also 
J.S. ROSEN [13]. 


61. There are 286 ninth-order equations of condition. However, this number 
of equations is greatly reduced if we make - very similar as in Parts I, II, 
and HI - the following assumptions: 
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<*14 =Q! 16 =1, Q!l 6 = 0, C x = C x =0, C 3 =C 2 =0, . . . , £7 =0? =0, 

A A A A A 

“ ^8 > Cq" °9 j •••> C 13 ~ C L3 j C 14 “ ^ > c 15 = C 16 “ 

3 3 1 = 341 =3q = • • » = 3 161 : 0, 052 = f^2 = 373 = . . . = 3lg2 = 0 

373 = 3q 3 = 3g3 = ... = 3l63 = 0, 3g4 = 3g4 = 3io4 = • • • = 3x©4 = 0 
and: 


(Al) 

Bn =}4 

= 2 , 3,4 

16 ) 

(A 2 ) 

T> 1 3 

Bvs - g 0i v 

(v = 2 , 3 , 4 , ..., 

16 ) 

(A 3 ) 

„ 1 4 

B13 - ^ 0l\) 

(v = 5 , 6 , 7 , ..., 

16 ) 

(A 4 ) 

ID 1 5 

p v* “ 5 a v 

(v = 8 , 9 , 10 , ... 

, 16 ) 


(137) 


(Al-v) 

(A 2 -V) 

(A 3 -v) 

(A 4 -v) 


(B) 



13 r 

C\)^V5 ”*"^14 ) 
Va 8 V 

3 146 ) 

3x55 + 3 165 / 

= c 5 ( 1 - 0 ( 5 ) = 0 

(B5) 

13 / 
^ Cv^V 6 + C 14 | 

3 146 ) 

■ 3l5S + 3l66 J 

= c%(l- 0 fe) = 0 

(B 6 ) 

13 

C\)0\77 + C^4 \ 
^*8 V ( 

3 147 ) 

3 137 + 3 367 / 

= O 7 ( 1 - 0(7 ) = 0 

(B7) 

13 , 

S C ^ + C X 4 ^ 
^s9 v l 

3 146 \ 

3l5B + 3l68 ' 

= Ca (1-Oe) 

(B 8 ) 

13 L 

X) c vBvje + c 14 < 

V=10 V 

3 146 ) 

3ig3 + 3le9 f 

= Cg (1-CKg ) 

(B9) 

13 

3 1410 ) 

3l510 + 3isi0 / 



^c v 3vio + | 

= c 10 (l~ttio) 

(BIO) 

13 

^ 2 c v 3 v u + c 14 { 

3l411 | 

3 1511 + 3 1611 J 

= c u (1-CKu) 

(Bll) 

C LB &L312 + c 14^ 

f 3 1412 ) 

(3 1512 + 3 1612 ; 

= Ci 2 (I - Q!l 2 ) 

(B12) 

C X 4< 

f 3l413 \ 

I3l513 + 3lsl3 J 

= ^13 (1 — 0fl3 ) 

(B13) 

c l4(3l514 + 01614) 

= C 14 (1-0(14 ) = 0 

(B14) 


C 14 01615 

= Ci4(l-0( 1B ) = C 14 

(B15) 
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C 13 

£ <^0^5 + Cl4 Pies = 0 (Cl) 

13 

53 C p,<^p / ^|JL6 + C 14 Pl66 = 0 (C2) 

13 

5Z C^CK^Bjj, 7 +^14 0X67 =0 (C3) 

S c ^v +cm ^ =o < c4) 

13 

(C) 53 C + Cx4 Pigg =0 (^°) 

13 

^Cp, 0^,0 JX7 +C14 0X67 =0 (C6) 

13 /M'l \ 13 

S C ^ a M»\S P ^ PvB / + Cl4 S 03 b v3 v s = 0 (C7) 

S c ^(s p p-V e ^) + S PxevPve = 0 (C8) 

i S c ^( i ' e ^^) +Ci4 = ° <c9) 


62. The assumptions of No. 61 reduce the equations of condition for our pair 
of eighth- and ninth-order Runge-Kutta formulas to the following two groups 
(D) and (E) in very much the same way as before in Parts I, II, or IE: 



(DO) 

(Dp) 
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SECTION XVII. A SOLUTION OF THE EQUATIONS OF 
CONDITION FOR THE RUNGE-KUTTA COEFFICIENTS 


63. We only briefly indicate how a solution of the equations of condition 
(Al), (A2), (A3), (A4), (B), (C), (D), and (E) can be obtained, since we 
follow a very similar procedure as in the earlier parts of our paper. 

The eight equations (Du) for the seven weight factors Crs , Cp , . . . , c 14 lead to a 
condition of compatibility for the coefficients a v . This condition and some 
other relations simplify considerably if we assume: 

ofe + c&f = 1, a 10 + an = 1, a x2 + a^ = 1 (138) 

Our condition of compatibility can then be written as: 

_ J_ . 1 . 42(5qfe -3) q 10 an a x3 «ia + 6 (4-7cife )Ja m a u + ai2a X3 ) - (5 - 9 oh ) , 139) 

8 3 2ob -1 70a X o an a X 2 an - 14(a-jo a-^ )+3 


If this relationship between a 8 , oq , a xo , a xx , a x2 , ffu holds, the weight factors 
C3 , Cp , . . . , c 14 can be determined from the first seven equations (Du). 


64. Besides (139) there are further restrictive conditions for the coefficients 
a v . Obviously the conditions (110), (111), (112) still hold in the case of our 
eighth-order formula. Furthermore, from equations (Al-8), (A2-8), (A3-8), 
(A4-8Xthe following condition results: 


1_ 20or a? ~ 15q R + 12qje 

q s - 5 a e 6 g£ a? -4a 0 +3c^ 


(140) 


From (110), (111), (112), (138), (139), (140), it follows that all a-coefficients 
can be expressed by q 4 , Og , 09 , a xo , a x2 . 

We shall show, however, in No. 63, that the coefficient afe is already deter- 
mined as soon as a e , a xo , ais are fixed. Therefore, q 4 , Og , q xo , a X 2 are 
the only free av-coefficients. 
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65. The 3-coefficients are computed in very much the same way as in Part 
ID. The expressions for Pa , P 33, P^, P 43, 053 , 054, % 3 , 064 , 0 6 s, 0?4 , 
075 , 076 are identical with those in Part in. 


Equations 


[ (A 1-8), (A2-8), (A3-8) 


085> 086 > 087 J 

) (A 1-9), (A2-9), (A3-9), (A4-9) | 

yield. 

095 > 096 y 097 > 0'98 1 

) (A1-10),(A2-10), (A3-10), (A4-10) 


0106 J 0107 > 01CS 9 0109 1 

\ (A1-11),(A2-11), (A3-11), (A4-11) 


0117 > 0118 j 0119 > 01110/ 


if we assume 


0105 — O, 0115 _ 0116 — 0 , 


(141) 


Assuming: 


0165 — 0145, 0166 — 0146, 0167 — 0147 


1 

(B5), (Cl), (C4) 


012 5 » 0135j 0145 

equations 

(B6), (C2), (C5) 

yield 

0126 9 0136 > 0148 


[ (B7), (C3), (C6) ] 


012? > 0137 j 0147 j 


From equations (Al-12), (A2-12), (A3-12), (A4-12), we can now obtain the 
coefficients P 13 ? , Pi$ , 0i2i O , 0isii . 


(142) 


66. Using the expressions for Pka as obtained from No. 65, the auxiliary 
quantities P 85 , P 95 , Piob , P 115 , P 125 are determined if ofe , a 8 , ofe , Oxo > c^is 
are given. 

From (B) it follows easily: 

C sP85 + °0 P 95 + CjjO Pl05 + Cxi Pxi5 + Cx2Pl25 + Cx 3 Pl35 +Cx4Pl46 

Equations (143), (El), (E3) can now be written in the form: 


73 



(144) 


C 13 Pl35 + c 14Pl45 S CvP^ 

k*8 

1 ^ 

c 13 0^13^X35 + c 14 Pi 46 = Cv^V^V 6 

, 12 

2 1 Y' 2 

Cl3 «X3Pl35 + C 14 P 145 =— ~ C v a v P V 5 

representing three equations for the two unknowns P 135 , P 146 . 

From (144) the following condition of compatibility can be obtained: 

12 1 
5Z c v( a i3 " a v)(l ” a v) ^vs = 3024 _ ^ q?13 (145) 

Equations (145) can be considered as a restrictive condition for ofe . Therefore, 
only a 4 , Oq, oi\o > ct\s can be considered as free parameters. 

67. The computation of the remaining 3 -coefficients is again straightforward 
and need hardly any explanation. 

From (A 1-13), (A2-13), (A3-13), (A4-13) , and the definition of P 135 - the 
quantity P 135 is already obtained from (144) - we obtain 3 L38 , 81 ®, 3 1310 , 81311 , 

3 1312 • 

The upper equations (B8) through (B13) yield 3 l4s , 3i® , P 1410 » Pi4u » 3 1412 , 

3 i 4 i 3 . This concludes the determination of the 3 -coefficients for our eighth - 
order formula. 



68 . We still have to determine the coefficients 615 land 3i6 \ for the ninth- 
order formula. 

From (E2) we can determine P 1wg . Since P^ =Pi 46 , equations (Al-16), 
(A2-16), (A3-16), (A4-16), the definition of P^s , and the definition of 
represent six equations for the six unknowns 8 * 33 , 81 ® , 8 i 6 io , 3i6ii> 81612 > 
81 B 13 , since the coefficients 8535 , Pies > 3i® are already determined by (142). 
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Equation (B14) suggests: 


Pl614 - Pl514 - 0 , 


(146) 


and equation (B15) yields: 


Pxas = 1 • (147) 

Finally, the coefficients Pi 5 \ can be determined by comparing the upper and 
the lower parts of equations (B) . 

All Runge-Kutta coefficients for our pair of eighth- and ninth-order formulas are 
now known if we still determine Cq , c 0 from equation (DO) and the coefficients 
Pyo (v = 1, 2, 3, . . . , 16) from the extended equations (F). If we submit the 
a-coefficients to the above restrictive conditions, the Runge-Kutta coefficients 
satisfy all our equations of condition. 
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SECTION XVIII. EXAMPLE FOR AN EIGHTH-ORDER 
RUNGE-KUTTA FORMULA 


69. In Table XII we present an example for an eighth-order Runge-Kutta 
formula RK8(9). Since cug , as obtained from equation (145), is in general 
no simple rational number anymore, we give in Table XII the decimal 
representation of the Runge-Kutta coefficients. The P- and c -coefficients 
not listed in Table XII are zero. The free parameters an , oie, ai 0 , 
were selected in such a way as to make the error coefficients of the leading 
error term of our Runge-Kutta formula as small as possible. 

The computation of the coefficients was performed on an IBM-7094 electronic 
computer in 40-digit arithmetic *. 


*The author is indebted to Mr. F.R. Calhoun, R-COMP-CSC, for 
developing the 40-digit package and making it available for these 
computations. 
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TABLE XII. RK 8 (9) 


O'! = 0.4436 8940 3764 9818 3109 5994 0428 1370 

Cfe = 0.6655 3410 5647 4727 4664 3991 0642 2055 

g? 3 = 0.9983 0115 8471 2091 1996 5986 5963 3083 

OM = 0.3155 0000 0000 0000 0000 0000 0000 0000 

or D = 0.5054 4100 9481 6906 8626 5161 2673 7384 

oi G = 0.1714 2857 1428 5714 2857 1428 5714 2857 

= 0.8285 7142 8571 4285 7142 8571 4285 7143 
eta = 0.6854 3966 1210 1156 2534 9537 6925 5586 

Os = 0.2487 8317 9680 6265 2069 7222 7456 0771 

a 10 = 0.1090 0000 0000 0000 0000 0000 0000 0000 

a 11 = 0.8910 0000 0000 0000 0000 0000 0000 0000 

ais = 0.3995 0000 0000 0000 0000 0000 0000 0000 

a 13 = 0.6005 0000 0000 0000 0000 0000 0000 0000 

<Xl4 = 1 

«15 = 0 

Ot 16 = 1 

8 10 = 0.4436 8940 3764 9818 3109 5994 0428 1370 

3so = 0.1663 8352 6411 8681 8666 0997 7660 5514 

Psi = 0.4991 5057 9235 6045 5998 2993 2981 6541 

830 = 0.2495 7528 9617 8022 7999 1496 6490 8271 

P32 = 0.7487 2586 8853 4068 3997 4489 9472 4812 

P40 = 0.2066 1891 1634 0060 2426 5567 1039 3185 

842 = 0.1770 7880 3779 8634 7040 3809 9728 8319 

=-0.6819 7715 4138 6949 4669 377 0 7681 5048 • 10 -1 
Pso = 0.1092 7823 1526 6640 8227 9038 9092 6157 

853 = 0.4021 5962 6423 6799 5421 9905 6369 0087 • 10" 2 

654 = 0.3921 4118 1690 7898 0444 3923 3017 4325 

8 60 = 0.9889 9281 4091 6466 5304 8447 6543 4355 • 10" 1 

8 g3 = 0.3513 8370 2279 6396 6951 2044 8735 6703 • 10 -2 

8 q4 = 0.1247 6099 9831 6001 6621 5206 2587 2489 

Pes =-0.5574 5546 8349 8979 9643 7429 0146 6348 • 10 -1 

P 70 =-0.3680 6865 2862 4220 3724 1531 0108 0691 

8 74 =-0.2227 3897 4694 7600 7645 0240 2094 4166 • 10 +1 

P ?5 = 0.1374 2908 2567 0291 0729 5656 9124 5744 • 10 +1 

8 7s = 0.2049 7390 0271 1160 3002 1593 5409 2206 • 10 +1 

8eo = 0.4546 7962 6413 4715 0077 3519 5060 3349 • 10 _1 
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TABLE XII. RK 8 (9) (Continued) 


P e6 = 0.3254 2131 
Pee = 0.2847 6660 
P^ = 0.9783 7801 
P 90 = 0.6084 2071 
Pg 5 -=-0.2118 4565 
Pse = 0.1959 6557 
p 97 =“0.4274 2640 
P 98 = 0.1743 4365 
Pjoo - 0.5405 9783 
P 1C6 = 0.1102 9325 
P 107 =-0.1256 5008 
P loe = 0.3679 0043 
Pi® -"0.5778 0542 
P no = 0.1273 2477 
P 117 = 0.1144 8305 
P U8 = 0.2877 3020 
P 119 = 0.5094 5379 
Pino ="0.1479 9682 
Pigo =-0.3652 6793 
P i;s = 0.8162 9896 
P 1S6 =-0.3860 7735 
Pis;, = 0.3086 2242 
Piss =-0.5807 7254 
P 1E0 = 0.3359 8659 
P 12l0 = 0.4106 6880 
8x211 s - 0.1184 0245 
Pi3o =-0.1237 5357 
Pjss =-0.2443 0768 
Piae = 0.5477 9568 
P 13? =-0.4441 3863 
Pias = 0.1001 3104 
Pus =-0.1499 5773 
Pi3io = 0.5894 6948 
01311 s 0.1738 0377 
01312= 0.2751 2330 
P 140 =-0.3526 0859 
0145 =-0.1839 6103 


7015 8914 7114 6774 
1385 2790 8888 1824 
6759 7915 2435 8683 
0626 2205 7051 0941 
7440 3700 7526 3252 
2661 7083 1957 4644 
3648 1760 3675 1448 
7368 1491 1965 3234 
2959 3191 7365 7857 
5978 2392 6539 2831 
5200 7255 6414 1477 
4775 8146 0136 3840 
7709 7207 3040 8406 
0636 6711 4646 6451 
0063 9310 5323 6588 
7096 9799 2776 2022 
4596 1136 3153 7358 
2443 7257 5900 2421 
8766 1674 0535 8485 
0123 1891 9777 8194 
6356 9350 6490 5176 
9246 0510 6450 4741 
5283 2060 2815 8293 
3288 8497 1493 1434 
4019 4995 8613 5496 
9723 5598 5520 6331 
9212 4514 3254 9790 
5513 5478 5358 7348 
9327 7865 6050 4365 
5334 1324 6374 9598 
8137 1326 6094 7926 
1020 5175 8447 1709 
5232 1701 3620 8245 
5034 2898 4877 6168 
6931 6673 0263 7586 
3883 3452 2700 5029 
1448 4827 0375 0441 


6964 8853 
2057 3687 
9727 1099 • 10 “ 3 
4520 5182 • 10 _1 
7525 1206 • 10 -1 
9066 2983 
3534 2899 • 10 -3 
5255 8189 • 10 " 1 
2411 1182 • 10 -1 
2764 8228 
6378 2250 • 10 “ 3 
4356 6339 • 10 -3 
2857 1866 • 10 _1 
8179 9160 
7572 1817 
0184 9198 
8507 9465 
4444 9640 
4439 4333 • 10 " 3 
2124 7030 • 10 -1 
9434 3215 
6602 5206 • 10 _1 
7473 3518 • 10 _1 
5136 2322 
2278 6417 
5615 4536 • 10 " 1 
9613 5669 • 10 +1 
6136 6763 • 10 +s 
2899 1173 
9656 9346 • 10 +1 
1785 1022 • 10 +s 
8507 3142 • 10 +3 
3965 1427 ■ 10 +1 
5744 0542 • 10 +1 
2286 0276 • 10 +3 
5887 5588 
9898 8231 
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TABLE XII. RK 8 (9) (Continued) 

0 146 =-0.6557 0189 4497 4164 5138 0068 7998 5251 

P 14 ? ="0.3908 6144 8804 3986 3435 0255 2024 1310 

0 x4e = 0.2679 4646 7128 5002 2936 5844 2327 1209 

3 14e =-0.1038 3022 9913 8249 0865 7698 5850 7427 • 10 +1 

3 i4io = 0.1667 2327 3242 5867 1664 7273 4616 8501 • 10 +1 

3i4U = 0.4955 1925 8553 1597 7067 7329 6707 1441 

3i4i2 = 0.1139 4001 1323 9706 3228 5867 3814 1784 • 10 +1 

8x413 = 0.5133 6696 4246 5861 3688 1990 9719 1534 • 10" 1 

3ieo = 0.1046 4847 3406 1481 0391 8730 0240 6755 • 10" 2 

Pies ="0.6716 3886 8449 9028 2237 7784 4617 8020 • 10 -2 

8ieb = 0.8182 8762 1894 2502 1265 3300 6524 8999 • 10~ 2 

Pisxo ="0.4264 0342 8644 8334 7277 1421 3808 7561 • 10 -2 

8i5ii = 0.2800 9029 4741 6893 6545 9763 3)15 3703 • lO^ 3 

Pibxb ="0.8783 5333 8762 3867 6639 0578 1314 5633 • 10~ 3 

01513 = 0.1025 4505 1108 2555 8084 2177 6966 4009 • 10" 1 

0160 ="0.1353 6550 7861 7406 7080 4421 6888 9966 • 10 +1 

0165 =“0.1839 6103 1448 4827 0375 0441 9898 8231 

Pies =“0.6557 0189 4497 4164 5138 0068 7998 5251 

0167 =-0.3908 6144 8804 3986 3435 0255 2024 1310 

P X6S = 0.2746 6285 5812 9992 5758 9622 0773 2989 

Pies =-0.1046 4851 7535 7191 5887 0351 8857 2676 • lO 4 ' 1 
Pisio = 0.1671 4967 6671 2315 5012 0044 8830 6588 • 10 4-1 

Pieii = 0.4952 3916 8258 4180 8131 1869 9074 0287 

01612 = 0.1148 1836 4662 7330 1905 2257 9595 4930 • 10 +1 

01613 = 0.4108 2191 3138 3305 5603 9813 2752 7525 • 10" 1 

01615 = 1 

e 0 = 0.3225 6083 5002 1624 9913 6129 0096 0247 • 10 _1 

c 8 = 0.2598 3725 2837 1540 3018 8870 2317 1963 

Cg = 0.9284 7805 9965 7702 7788 0637 1430 2190 • 10“ 

Cio = 0.1645 2339 5147 6434 2891 6477 3184 2800 

Cu = 0.1766 5951 6378 6007 4367 0842 9839 7547 

c L2 = 0.2392 0102 3203 5275 9374 1089 3332 0941 

Ci3 = 0.3948 4274 6042 0285 3746 7521 1882 9325 • 10" 2 

Ci 4 = 0.3072 6495 4758 6064 0406 3683 0552 2124 • 10 _1 

TE = Ci 4 (f 0 + fi4 - fis - fi 6 ) h (148) 
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SECTION XIX. NUMERICAL COMPARISON WITH OTHER 
EIGHTH-ORDER RUNGE-KUTTA FORMULAS 


70. We compare our formula RK8(9) with SHANKS's eighth-order formula 
(8-12) ([13], p.34) which is the only other eighth-order Runge-Kutta formula 
known to us . We apply both formulas to our problem (53) in exactly the same 
way as in Parts I, II, and III using RICHARDSON'S principle as stepsize 
control procedure for SHANKS's formula. Table XIII shows our results. 


TABLE XIII. COMPARISON OF EIGHTH-ORDER METHODS FOR EXAMPLE (53) 



Number of 

Results for x = 5 and Tolerance 10 v 

A „ 

Method 

Substitutions 

^Number of 

Total Number 
of Evaluations 

Running Time 
on rBM-7094 
(min) 

Accumulated Errors in y and z ’ 


per Step 

Steps 

Ay 

£ z 

SHANKS 

23 

G94 

15 962 

1.65 

-0.1710. 10" ir ' 

-0.4G4G • 10 _V5 

RK5(9) 

17 

510 

S 670 

0.91 

-0.177G.10 -14 

-0.3553- 10"^ 


The table shows that, in this example, the computer running time for our 
formula RK8(9) is about 55% of the running time of SHANKS's eighth-order 
formula . 
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